Changeset 828


Ignore:
Timestamp:
Feb 21, 2012 12:00:36 PM (13 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:
10 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r826 r828  
    44# Current revisions:
    55# ------------------
    6 #
     6# init_particles depends on lpm_collision_kernels
    77#
    88# Former revisions:
     
    6666
    6767RCS = advec_particles.f90 advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \
    68         advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \
    69         advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \
    70         advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \
    71         average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \
    72         calc_liquid_water_content.f90 calc_precipitation.f90 \
    73         calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \
    74         check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \
    75         coriolis.f90 cpu_log.f90 cpu_statistics.f90 data_log.f90 \
    76         data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \
    77         data_output_ptseries.f90 data_output_spectra.f90 data_output_tseries.f90 \
    78         data_output_2d.f90 data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 \
    79         diffusion_u.f90 diffusion_v.f90 diffusion_w.f90 diffusivities.f90 \
    80         disturb_field.f90 disturb_heatflux.f90 eqn_state_seawater.f90 \
    81         exchange_horiz.f90 exchange_horiz_2d.f90 \
    82         fft_xy.f90 flow_statistics.f90 global_min_max.f90 \
    83         header.f90 impact_of_latent_heat.f90 inflow_turbulence.f90 \
    84         init_1d_model.f90 init_3d_model.f90 init_advec.f90 \
    85         init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 init_grid.f90 \
    86         init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \
    87         init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
     68        advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \
     69        advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \
     70        advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \
     71        average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \
     72        calc_liquid_water_content.f90 calc_precipitation.f90 \
     73        calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \
     74        check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \
     75        coriolis.f90 cpu_log.f90 cpu_statistics.f90 data_log.f90 \
     76        data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \
     77        data_output_ptseries.f90 data_output_spectra.f90 \
     78        data_output_tseries.f90 data_output_2d.f90 data_output_3d.f90 \
     79        diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 diffusion_v.f90 \
     80        diffusion_w.f90 diffusivities.f90 disturb_field.f90 \
     81        disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \
     82        exchange_horiz_2d.f90 fft_xy.f90 flow_statistics.f90 \
     83        global_min_max.f90 header.f90 impact_of_latent_heat.f90 \
     84        inflow_turbulence.f90 init_1d_model.f90 init_3d_model.f90 \
     85        init_advec.f90 init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 \
     86        init_grid.f90 init_masks.f90 init_ocean.f90 init_particles.f90 \
     87        init_pegrid.f90 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
    8888        interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \
    89         local_stop.f90 local_system.f90 local_tremain.f90 local_tremain_ini.f90 \
    90         lpm_collision_kernels.f90 message.f90 modules.f90 netcdf.f90 \
    91         package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \
    92         plant_canopy_model.f90 poisfft.f90 \
     89        local_stop.f90 local_system.f90 local_tremain.f90 \
     90        local_tremain_ini.f90 lpm_collision_kernels.f90 message.f90 \
     91        modules.f90 netcdf.f90 package_parin.f90 palm.f90 parin.f90 \
     92        particle_boundary_conds.f90 plant_canopy_model.f90 poisfft.f90 \
    9393        poisfft_hybrid.f90 poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \
    9494        production_e.f90 prognostic_equations.f90 random_function.f90 \
    9595        random_gauss.f90 read_3d_binary.f90 read_var_list.f90 run_control.f90 \
    9696        set_particle_attributes.f90 set_slicer_attributes_dvrp.f90 \
    97         singleton.f90 sor.f90 spline_x.f90 \
     97        singleton.f90 sor.f90 spline_x.f90 \
    9898        spline_y.f90 spline_z.f90 subsidence.f90 \
    99         sum_up_3d_data.f90 surface_coupler.f90 \
     99        sum_up_3d_data.f90 surface_coupler.f90 \
    100100        swap_timelevel.f90 temperton_fft.f90 time_integration.f90 \
    101         time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \
    102         transpose.f90 user_3d_data_averaging.f90 user_actions.f90 \
     101        time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \
     102        transpose.f90 user_3d_data_averaging.f90 user_actions.f90 \
    103103        user_additional_routines.f90 user_advec_particles.f90 \
    104104        user_check_data_output.f90 user_check_data_output_pr.f90 \
    105105        user_check_parameters.f90 user_data_output_2d.f90 \
    106         user_data_output_3d.f90 user_data_output_dvrp.f90 user_data_output_mask.f90 \
    107         user_define_netcdf_grid.f90 user_dvrp_coltab.f90 user_header.f90 \
    108         user_init.f90 user_init_3d_model.f90 user_init_grid.f90 \
    109         user_init_particles.f90 user_init_plant_canopy.f90 \
    110         user_last_actions.f90 user_module.f90 user_parin.f90 \
    111         user_particle_attributes.f90 user_read_restart_data.f90 \
     106        user_data_output_3d.f90 user_data_output_dvrp.f90 \
     107        user_data_output_mask.f90 user_define_netcdf_grid.f90 \
     108        user_dvrp_coltab.f90 user_header.f90 user_init.f90 \
     109        user_init_3d_model.f90 user_init_grid.f90 user_init_particles.f90 \
     110        user_init_plant_canopy.f90 user_last_actions.f90 user_module.f90 \
     111        user_parin.f90 user_particle_attributes.f90 user_read_restart_data.f90 \
    112112        user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
    113         write_3d_binary.f90 write_compressed.f90 \
    114         write_var_list.f90
     113        write_3d_binary.f90 write_compressed.f90 write_var_list.f90
    115114
    116115OBJS =  advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \
     
    252251init_masks.o: modules.o
    253252init_ocean.o: modules.o eqn_state_seawater.o
    254 init_particles.o: modules.o random_function.o
     253init_particles.o: modules.o lpm_collision_kernels.o random_function.o
    255254init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o
    256255init_pt_anomaly.o: modules.o
  • palm/trunk/SOURCE/advec_particles.f90

    r826 r828  
    44! Current revisions:
    55! ------------------
    6 !
     6! fast hall/wang kernels with fixed radius/dissipation classes added,
     7! particle feature color renamed class, routine colker renamed
     8! recalculate_kernel,
     9! lower limit for droplet radius changed from 1E-7 to 1E-8
     10!
     11! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
    712!
    813! Former revisions:
     
    144149    IMPLICIT NONE
    145150
    146     INTEGER ::  agp, deleted_particles, deleted_tails, i, ie, ii, inc,        &
    147                 internal_timestep_count, is, j, &
    148                 jj, js, jtry, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi,   &
     151    INTEGER ::  agp, deleted_particles, deleted_tails, eclass, i, ie, ii, inc, &
     152                internal_timestep_count, is, j, jj, js, jtry, k, kk, kw, m, n, &
     153                nc, nd, nn, num_gp, pse, psi, rclass_l, rclass_s,  &
    149154                tlength, trlp_count, trlp_count_sum, trlp_count_recv,          &
    150155                trlp_count_recv_sum, trlpt_count, trlpt_count_recv,            &
     
    173178                lagr_timescale, lw_max, mean_r, new_r, p_int, pt_int, &
    174179                pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, &
    175                 ql_int_u, random_gauss, r_ros, r_ros_ini, sigma, sl_r3, sl_r4, &
    176                 t_int, u_int, u_int_l, u_int_u,vv_int, v_int, v_int_l, &
    177                 v_int_u, w_int, w_int_l, w_int_u, x, y
     180                ql_int_u, random_gauss, r_ros, r_ros_ini, &
     181                sigma, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,vv_int, &
     182                v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u, x, y
    178183!
    179184!-- Parameters for Rosenbrock method
     
    190195
    191196    REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
    192     REAL, DIMENSION(:,:), ALLOCATABLE ::  kern
    193197
    194198    REAL    ::  location(1:30,1:3)
     
    356360                             thermal_conductivity_l +                         &
    357361                             rho_l * r_v * t_int / diff_coeff_l / e_s )
    358              IF ( arg < 1.0E-14 )  THEN
    359                 new_r = 1.0E-7
     362             IF ( arg < 1.0E-16 )  THEN
     363                new_r = 1.0E-8
    360364             ELSE
    361365                new_r = SQRT( arg )
     
    550554             CALL message( 'advec_particles', 'PA0144', 2, 2, -1, 6, 1 )
    551555          ENDIF
    552           particles(n)%radius = new_r
    553556
    554557!
    555558!--       Sum up the total volume of liquid water (needed below for
    556559!--       re-calculating the weighting factors)
    557           ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor *             &
    558                                       particles(n)%radius**3
     560          ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r**3
     561
     562          particles(n)%radius = new_r
     563
     564!
     565!--       Determine radius class of the particle needed for collision
     566          IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  use_kernel_tables ) &
     567          THEN
     568             particles(n)%class = ( LOG( new_r ) - rclass_lbound ) /  &
     569                                  ( rclass_ubound - rclass_lbound ) * &
     570                                  radius_classes
     571             particles(n)%class = MIN( particles(n)%class, radius_classes )
     572             particles(n)%class = MAX( particles(n)%class, 1 )
     573          ENDIF
     574
    559575       ENDDO
    560576       CALL cpu_log( log_point_s(42), 'advec_part_cond', 'stop' )
     
    601617                   pse = psi + prt_count(k,j,i)-1
    602618
    603                    IF ( hall_kernel  .OR.  wang_kernel )  THEN
    604 
    605                       ALLOCATE( kern(psi:pse,psi:pse) )
    606 
    607 !
    608 !--                   Calculate collision kernel for all particles in grid box
    609                       CALL colker( i, j, k, kern )
    610 !
    611 !--                   Collison kernel is calculated in cm**3/s but needed
    612 !--                   in m**3/s
    613                       kern = kern * 1.0E-6   ! to be moved to colker
     619!
     620!--                Now apply the different kernels
     621                   IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
     622                        use_kernel_tables )  THEN
     623!
     624!--                   Fast method with pre-calculated efficiencies for
     625!--                   discrete radius- and dissipation-classes.
     626!
     627!--                   Determine dissipation class index of this gridbox
     628                      eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * &
     629                                    dissipation_classes ) + 1
     630                      IF ( hall_kernel .OR. diss(k,j,i) * 1.0E4 < 0.001 )  THEN
     631                         eclass = 0   ! Hall kernel is used
     632                      ELSE
     633                         eclass = MIN( dissipation_classes, eclass )
     634                      ENDIF
    614635
    615636                      DO  n = pse, psi+1, -1
     
    617638                         integral = 0.0
    618639                         lw_max   = 0.0
    619 
     640                         rclass_l = particles(n)%class
    620641!
    621642!--                      Calculate growth of collector particle radius using all
     
    623644                         DO  is = psi, n-1
    624645
    625                             integral = integral +                              &
    626                                        ( particles(is)%radius**3 * kern(n,is) *&
     646                            rclass_s = particles(is)%class
     647                            integral = integral +                             &
     648                                       ( particles(is)%radius**3 *            &
     649                                         ckernel(rclass_l,rclass_s,eclass) *  &
    627650                                         particles(is)%weight_factor )
    628651!
     
    653676
    654677                         ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
    655 
     678!-- can this case really happen??
    656679                            DO  is = psi, n-1
    657680
     
    667690                            DO  is = psi, n-1
    668691
     692                               rclass_s = particles(is)%class
    669693                               particles(is)%weight_factor = &
    670694                                  particles(is)%weight_factor - &
    671                                   ( ( kern(n,is) * particles(is)%weight_factor &
     695                                  ( ( ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor &
    672696                                    / integral ) * delta_v )
    673697
     
    699723                      ENDDO
    700724
    701                       DEALLOCATE( kern )
     725                   ELSEIF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
     726                            .NOT. use_kernel_tables )  THEN
     727!
     728!--                   Collision efficiencies are calculated for every new
     729!--                   grid box. First, allocate memory for kernel table.
     730!--                   Third dimension is 1, because table is re-calculated for
     731!--                   every new dissipation value.
     732                      ALLOCATE( ckernel(prt_start_index(k,j,i):                &
     733                               prt_start_index(k,j,i)+prt_count(k,j,i)-1,      &
     734                               prt_start_index(k,j,i):                         &
     735                               prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) )
     736!
     737!--                   Now calculate collision efficiencies for this box
     738                      CALL recalculate_kernel( i, j, k )
     739
     740                      DO  n = pse, psi+1, -1
     741
     742                         integral = 0.0
     743                         lw_max   = 0.0
     744!
     745!--                      Calculate growth of collector particle radius using all
     746!--                      droplets smaller than current droplet
     747                         DO  is = psi, n-1
     748
     749                            integral = integral + particles(is)%radius**3 * &
     750                                                  ckernel(n,is,1) *         &
     751                                                  particles(is)%weight_factor
     752!
     753!--                         Calculate volume of liquid water of the collected
     754!--                         droplets which is the maximum liquid water available
     755!--                         for droplet growth   
     756                            lw_max = lw_max + ( particles(is)%radius**3 * &
     757                                                particles(is)%weight_factor )
     758
     759                         ENDDO
     760
     761!
     762!--                      Change in radius of collector droplet due to collision
     763                         delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
     764                                   integral * dt_3d * ddx * ddy / dz
     765
     766!
     767!--                      Change in volume of collector droplet due to collision
     768                         delta_v = particles(n)%weight_factor                  &
     769                                   * ( ( particles(n)%radius + delta_r )**3    &
     770                                       - particles(n)%radius**3 )
     771
     772                         IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
     773!-- replace by message call
     774                            PRINT*, 'Particle has grown to much because the',  &
     775                                    ' change of volume of particle is larger', &
     776                                    ' than liquid water available!'
     777
     778                         ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
     779!-- can this case really happen??
     780                            DO  is = psi, n-1
     781
     782                               particles(is)%weight_factor = 0.0
     783                               particle_mask(is)  = .FALSE.
     784                               deleted_particles = deleted_particles + 1
     785
     786                            ENDDO
     787
     788                         ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
     789!
     790!--                         Calculate new weighting factor of collected droplets
     791                            DO  is = psi, n-1
     792
     793                               particles(is)%weight_factor =                   &
     794                                      particles(is)%weight_factor -            &
     795                                      ( ckernel(n,is,1) / integral * delta_v * &
     796                                        particles(is)%weight_factor )
     797
     798                               IF ( particles(is)%weight_factor < 0.0 )  THEN
     799                                     WRITE( message_string, * ) 'negative ',   &
     800                                                        'weighting factor: ',  &
     801                                                     particles(is)%weight_factor
     802                                     CALL message( 'advec_particles', '', 2, 2,&
     803                                                   -1, 6, 1 )
     804
     805                               ELSEIF ( particles(is)%weight_factor == 0.0 )  &
     806                               THEN
     807
     808                                   particles(is)%weight_factor = 0.0
     809                                   particle_mask(is)  = .FALSE.
     810                                   deleted_particles = deleted_particles + 1
     811
     812                               ENDIF
     813
     814                            ENDDO
     815
     816                         ENDIF
     817
     818                         particles(n)%radius = particles(n)%radius + delta_r
     819                         ql_vp(k,j,i) = ql_vp(k,j,i) + &
     820                                                    particles(n)%weight_factor &
     821                                                    * particles(n)%radius**3
     822
     823                      ENDDO
     824
     825                      DEALLOCATE( ckernel )
    702826
    703827                   ELSEIF ( palm_kernel )  THEN
     
    23762500                      particle_tail_coordinates(1,2,nn) = particles(n)%y
    23772501                      particle_tail_coordinates(1,3,nn) = particles(n)%z
    2378                       particle_tail_coordinates(1,4,nn) = particles(n)%color
     2502                      particle_tail_coordinates(1,4,nn) = particles(n)%class
    23792503                      particles(n)%tailpoints = 1
    23802504                      IF ( minimum_tailpoint_distance /= 0.0 )  THEN
     
    23822506                         particle_tail_coordinates(2,2,nn) = particles(n)%y
    23832507                         particle_tail_coordinates(2,3,nn) = particles(n)%z
    2384                          particle_tail_coordinates(2,4,nn) = particles(n)%color
     2508                         particle_tail_coordinates(2,4,nn) = particles(n)%class
    23852509                         particle_tail_coordinates(1:2,5,nn) = 0.0
    23862510                         particles(n)%tailpoints = 2
     
    37433867
    37443868!
    3745 !-- Set particle attributes
    3746     CALL set_particle_attributes
     3869!-- Set particle attributes.
     3870!-- Feature is not available if collision is activated, because the respective
     3871!-- particle attribute (class) is then used for storing the particle radius
     3872!-- class.
     3873    IF ( collision_kernel == 'none' )  CALL set_particle_attributes
    37473874
    37483875!
    37493876!-- Set particle attributes defined by the user
    37503877    CALL user_particle_attributes
    3751 !    WRITE ( 9, * ) '*** advec_particles: ##10'
    3752 !    CALL local_flush( 9 )
    3753 !    DO  n = 1, number_of_particles
    3754 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    3755 !       THEN
    3756 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3757 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    3758 !          CALL local_flush( 9 )
    3759 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3760 !       ENDIF
    3761 !    ENDDO
    37623878
    37633879!
     
    38103926             particle_tail_coordinates(1,2,nn) = particles(n)%y
    38113927             particle_tail_coordinates(1,3,nn) = particles(n)%z
    3812              particle_tail_coordinates(1,4,nn) = particles(n)%color
     3928             particle_tail_coordinates(1,4,nn) = particles(n)%class
    38133929!             WRITE ( 9, * ) '*** advec_particles: ##10.4'
    38143930!             CALL local_flush( 9 )
     
    40584174    CALL handle_netcdf_error( 'output_particles_netcdf', 15 )
    40594175
    4060     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%color,      &
     4176    nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%class,      &
    40614177                            start = (/ 1, prt_time_count /),                  &
    40624178                            count = (/ maximum_number_of_particles /) )
  • palm/trunk/SOURCE/check_parameters.f90

    r826 r828  
    44! Current revisions:
    55! -----------------
    6 !
     6! check of collision_kernel extended
    77!
    88! Former revisions:
     
    716716    SELECT CASE ( TRIM( collision_kernel ) )
    717717
    718        CASE ( 'hall' )
     718       CASE ( 'hall', 'hall_fast' )
    719719          hall_kernel = .TRUE.
    720720
     
    722722          palm_kernel = .TRUE.
    723723
    724        CASE ( 'wang' )
     724       CASE ( 'wang', 'wang_fast' )
    725725          wang_kernel = .TRUE.
    726726
     
    734734
    735735    END SELECT
     736    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
    736737
    737738
  • palm/trunk/SOURCE/data_output_dvrp.f90

    r392 r828  
    3232! Current revisions:
    3333! -----------------
    34 !
     34! particle feature color renamed class
    3535!
    3636! Former revisions:
     
    226226                   p_y(k)   = particles(n)%y * superelevation_y
    227227                   p_z(k)   = particles(n)%z * superelevation
    228                    p_c(k)   = particles(n)%color
     228                   p_c(k)   = particles(n)%class
    229229                ENDIF
    230230             ENDDO
  • palm/trunk/SOURCE/header.f90

    r826 r828  
    44! Current revisions:
    55! -----------------
    6 !
     6! output of dissipation_classes + radius_classes
    77!
    88! Former revisions:
     
    14831483       IF ( collision_kernel /= 'none' )  THEN
    14841484          WRITE ( io, 435 )  TRIM( collision_kernel )
    1485        ELSE
    1486           WRITE ( io, 436 )
     1485          IF ( collision_kernel(6:9) == 'fast' )  THEN
     1486             WRITE ( io, 436 )  radius_classes, dissipation_classes
     1487          ENDIF
     1488       ELSE
     1489          WRITE ( io, 437 )
    14871490       ENDIF
    14881491    ENDIF
     
    19641967                 ' droplets < 1.0E-6 m')
    19651968435 FORMAT ('    Droplet collision is handled by ',A,'-kernel')
    1966 436 FORMAT ('    Droplet collision is switched off')
     1969436 FORMAT ('       Fast kernel with fixed radius- and dissipation classes ', &
     1970                    'are used'/ &
     1971            '          number of radius classes:       ',I3,'    interval ', &
     1972                       '[1.0E-6,2.0E-4] m'/ &
     1973            '          number of dissipation classes:   ',I2,'    interval ', &
     1974                       '[0,1000] cm**2/s**3')
     1975437 FORMAT ('    Droplet collision is switched off')
    19671976450 FORMAT (//' LES / Turbulence quantities:'/ &
    19681977              ' ---------------------------'/)
  • palm/trunk/SOURCE/init_particles.f90

    r826 r828  
    44! Current revisions:
    55! -----------------
    6 !
     6! call of init_kernels, particle feature color renamed class
    77!
    88! Former revisions:
     
    7272    USE grid_variables
    7373    USE indices
     74    USE lpm_collision_kernels_mod
    7475    USE particle_attributes
    7576    USE pegrid
     
    148149       IF ( pdz(j) == 9999999.9  .OR.  pdz(j) == 0.0 )  pdz(j) = pdz(j-1)
    149150    ENDDO
     151
     152!
     153!-- Initialize collision kernels
     154    IF ( collision_kernel /= 'none' )  CALL init_kernels
    150155
    151156!
     
    343348                            particles(n)%radius      = particle_groups(i)%radius
    344349                            particles(n)%weight_factor =initial_weighting_factor
    345                             particles(n)%color         = 1
     350                            particles(n)%class         = 1
    346351                            particles(n)%group         = i
    347352                            particles(n)%tailpoints    = 0
     
    470475             IF ( MINVAL( particles(1:n)%dvrp_psize  ) ==     &
    471476                  MAXVAL( particles(1:n)%dvrp_psize  )  .AND. &
    472                   MINVAL( particles(1:n)%color ) ==     &
    473                   MAXVAL( particles(1:n)%color ) )  THEN
     477                  MINVAL( particles(1:n)%class ) ==     &
     478                  MAXVAL( particles(1:n)%class ) )  THEN
    474479                uniform_particles_l = .TRUE.
    475480             ELSE
     
    532537                particle_tail_coordinates(1,2,nn) = particles(n)%y
    533538                particle_tail_coordinates(1,3,nn) = particles(n)%z
    534                 particle_tail_coordinates(1,4,nn) = particles(n)%color
     539                particle_tail_coordinates(1,4,nn) = particles(n)%class
    535540                particles(n)%tailpoints = 1
    536541                IF ( minimum_tailpoint_distance /= 0.0 )  THEN
     
    538543                   particle_tail_coordinates(2,2,nn) = particles(n)%y
    539544                   particle_tail_coordinates(2,3,nn) = particles(n)%z
    540                    particle_tail_coordinates(2,4,nn) = particles(n)%color
     545                   particle_tail_coordinates(2,4,nn) = particles(n)%class
    541546                   particle_tail_coordinates(1:2,5,nn) = 0.0
    542547                   particles(n)%tailpoints = 2
  • palm/trunk/SOURCE/lpm_collision_kernels.f90

    r826 r828  
    1 MODULE lpm_collision_kernels_mod
     1 MODULE lpm_collision_kernels_mod
    22
    33!------------------------------------------------------------------------------!
    44! Current revisions:
    55! -----------------
    6 !
     6! code has been completely reformatted, routine colker renamed
     7! recalculate_kernel,
     8! routine init_kernels added, radius is now communicated to the collision
     9! routines by array radclass
     10!
     11! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
    712!
    813! Former revisions:
     
    2429! Description:
    2530! ------------
    26 ! This routine calculates the effect of (SGS) turbulence on the collision
    27 ! efficiency of droplets.
    28 ! It is based on the original kernel developed by Wang (...)
     31! This module calculates collision efficiencies either due to pure gravitational
     32! effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
     33! including the effects of (SGS) turbulence (Wang kernel, see Wang and
     34! Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been
     35! provided by L.-P. Wang but is substantially reformatted and speed optimized
     36! here.
     37!
     38! ATTENTION:
     39! Physical quantities (like g, densities, etc.) used in this module still
     40! have to be adjusted to those values used in the main PALM code.
     41! Also, quantities in CGS-units should be converted to SI-units eventually.
    2942!------------------------------------------------------------------------------!
    3043
     
    3346    USE constants
    3447    USE particle_attributes
     48    USE pegrid
     49
    3550
    3651    IMPLICIT NONE
     
    3853    PRIVATE
    3954
    40     PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi
    41 
    42     INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
    43 
    44     REAL, SAVE :: epsilon, eps2, urms, urms2
    45 
    46     REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
    47     REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
     55    PUBLIC  ckernel, init_kernels,  rclass_lbound, rclass_ubound, &
     56            recalculate_kernel           
     57
     58    REAL ::  epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2
     59
     60    REAL, DIMENSION(:),   ALLOCATABLE ::  epsclass, radclass, winf
     61    REAL, DIMENSION(:,:), ALLOCATABLE ::  ec, ecf, gck, hkernel, hwratio
     62    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  ckernel
     63
     64    SAVE
    4865
    4966!
    5067!-- Public interfaces
    51     INTERFACE colker
    52        MODULE PROCEDURE colker
    53     END INTERFACE colker
    54 
    55     INTERFACE effic
    56        MODULE PROCEDURE effic
    57     END INTERFACE effic
    58 
    59     INTERFACE fallg
    60        MODULE PROCEDURE fallg
    61     END INTERFACE fallg
    62 
    63     INTERFACE phi
    64        MODULE PROCEDURE phi
    65     END INTERFACE phi
    66 
    67     INTERFACE turbsd
    68        MODULE PROCEDURE turbsd
    69     END INTERFACE turbsd
    70 
    71     INTERFACE turb_enhance_eff
    72        MODULE PROCEDURE turb_enhance_eff
    73     END INTERFACE turb_enhance_eff
    74 
    75     INTERFACE zhi
    76        MODULE PROCEDURE zhi
    77     END INTERFACE zhi
     68    INTERFACE init_kernels
     69       MODULE PROCEDURE init_kernels
     70    END INTERFACE init_kernels
     71
     72    INTERFACE recalculate_kernel
     73       MODULE PROCEDURE recalculate_kernel
     74    END INTERFACE recalculate_kernel
     75
    7876
    7977    CONTAINS
    8078
    81 !------------------------------------------------------------------------------!
    82 ! SUBROUTINE for calculation of collision kernel
    83 !------------------------------------------------------------------------------!
    84     SUBROUTINE colker( i1, j1, k1, kernel )
     79
     80    SUBROUTINE init_kernels
     81!------------------------------------------------------------------------------!
     82! Initialization of the collision efficiency matrix with fixed radius and
     83! dissipation classes, calculated at simulation start only.
     84!------------------------------------------------------------------------------!
     85
     86       IMPLICIT NONE
     87
     88       INTEGER ::  i, j, k
     89
     90
     91!
     92!--    Calculate collision efficiencies for fixed radius- and dissipation
     93!--    classes
     94       IF ( collision_kernel(6:9) == 'fast' )  THEN
     95
     96          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,               &
     97                    0:dissipation_classes), epsclass(1:dissipation_classes), &
     98                    radclass(1:radius_classes) )
     99
     100!
     101!--       Calculate the radius class bounds with logarithmic distances
     102!--       in the interval [1.0E-6, 2.0E-4] m
     103          rclass_lbound = LOG( 1.0E-6 )
     104          rclass_ubound = LOG( 2.0E-4 )
     105          radclass(1)   = 1.0E-6
     106          DO  i = 2, radius_classes
     107             radclass(i) = EXP( rclass_lbound +                                &
     108                                ( rclass_ubound - rclass_lbound ) * ( i-1.0 ) /&
     109                                ( radius_classes - 1.0 ) )
     110!             IF ( myid == 0 )  THEN
     111!                PRINT*, 'i=', i, ' r = ', radclass(i)*1.0E6
     112!             ENDIF
     113          ENDDO
     114!
     115!--       Collision routines expect radius to be in cm
     116          radclass = radclass * 100.0
     117
     118!
     119!--       Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3
     120          DO  i = 1, dissipation_classes
     121             epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes
     122!             IF ( myid == 0 )  THEN
     123!                PRINT*, 'i=', i, ' eps = ', epsclass(i)
     124!             ENDIF
     125          ENDDO
     126!
     127!--       Calculate collision efficiencies of the Wang/ayala kernel
     128          ALLOCATE( ec(1:radius_classes,1:radius_classes),  &
     129                    ecf(1:radius_classes,1:radius_classes), &
     130                    gck(1:radius_classes,1:radius_classes), &
     131                    winf(1:radius_classes) )
     132
     133          DO  k = 1, dissipation_classes
     134
     135             epsilon = epsclass(k)
     136             urms    = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 )
     137
     138             CALL turbsd
     139             CALL turb_enhance_eff
     140             CALL effic
     141
     142             DO  j = 1, radius_classes
     143                DO  i = 1, radius_classes
     144                   ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j)
     145                ENDDO
     146             ENDDO
     147
     148          ENDDO
     149
     150!
     151!--       Calculate collision efficiencies of the Hall kernel
     152          ALLOCATE( hkernel(1:radius_classes,1:radius_classes), &
     153                    hwratio(1:radius_classes,1:radius_classes) )
     154
     155          CALL fallg
     156          CALL effic
     157
     158          DO  j = 1, radius_classes
     159             DO  i =  1, radius_classes
     160                hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 &
     161                                  * ec(i,j) * ABS( winf(j) - winf(i) )
     162                ckernel(i,j,0) = hkernel(i,j)  ! hall kernel stored on index 0
     163              ENDDO
     164          ENDDO
     165
     166!
     167!--       Test output of efficiencies
     168          IF ( j == -1 )  THEN
     169
     170             PRINT*, '*** Hall kernel'
     171             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
     172             DO  j = 1, radius_classes
     173                WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes )
     174             ENDDO
     175
     176             DO  k = 1, dissipation_classes
     177                DO  i = 1, radius_classes
     178                   DO  j = 1, radius_classes
     179                      IF ( hkernel(i,j) == 0.0 )  THEN
     180                         hwratio(i,j) = 9999999.9
     181                      ELSE
     182                         hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j)
     183                      ENDIF
     184                   ENDDO
     185                ENDDO
     186
     187                PRINT*, '*** epsilon = ', epsclass(k)
     188                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
     189                DO  j = 1, radius_classes
     190!                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes )
     191                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes )
     192                ENDDO
     193             ENDDO
     194
     195          ENDIF
     196
     197          DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
     198
     199          ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
     200
     201       ELSEIF( collision_kernel == 'hall'  .OR.  collision_kernel == 'wang' ) &
     202       THEN
     203!
     204!--       Initial settings for Hall- and Wang-Kernel
     205!--       To be done: move here parts from turbsd, fallg, ecoll, etc.
     206       ENDIF
     207
     208    END SUBROUTINE init_kernels
     209
     210
     211!------------------------------------------------------------------------------!
     212! Calculation of collision kernels during each timestep and for each grid box
     213!------------------------------------------------------------------------------!
     214    SUBROUTINE recalculate_kernel( i1, j1, k1 )
    85215
    86216       USE arrays_3d
     
    94224       IMPLICIT NONE
    95225
    96        INTEGER ::  i, i1, j, j1, k1
    97 
    98        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
    99 
    100 !       CALL cpu_log( log_point_s(46), 'colker', 'start' )
    101 
    102        ip = i1
    103        jp = j1
    104        kp = k1
    105 
    106        pstart = prt_start_index(kp,jp,ip)
    107        pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
    108        psum   = prt_count(kp,jp,ip)
    109 
    110        ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
    111 
    112        IF ( wang_kernel )  THEN
    113 
    114           ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
    115 
    116           epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
    117           urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
    118 
    119           IF ( epsilon <= 0.001 )  THEN
    120 
    121              CALL fallg
    122              CALL effic
    123 
    124              DO  j = pstart, pend
    125                 DO  i =  pstart, pend
    126                    kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
    127                                         particles(i)%radius * 100.0 )**2 &
    128                                     * ec(i,j) * ABS( winf(j) - winf(i) )
    129                 ENDDO
     226       INTEGER ::  i, i1, j, j1, k1, pend, pstart
     227
     228
     229       pstart = prt_start_index(k1,j1,i1)
     230       pend   = prt_start_index(k1,j1,i1) + prt_count(k1,j1,i1) - 1
     231       radius_classes = prt_count(k1,j1,i1)
     232
     233       ALLOCATE( ec(1:radius_classes,1:radius_classes), &
     234                 radclass(1:radius_classes), winf(1:radius_classes) )
     235
     236!
     237!--    Store particle radii on the radclass array. Collision routines
     238!--    expect radii to be in cm.
     239       radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0
     240
     241       epsilon = diss(k1,j1,i1) * 1.0E4   ! dissipation rate in cm**2/s**-3
     242       urms    = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 )
     243
     244       IF ( wang_kernel  .AND.  epsilon > 0.001 )  THEN
     245!
     246!--       Call routines to calculate efficiencies for the Wang kernel
     247          ALLOCATE( gck(1:radius_classes,1:radius_classes), &
     248                    ecf(1:radius_classes,1:radius_classes) )
     249
     250          CALL turbsd
     251          CALL turb_enhance_eff
     252          CALL effic
     253
     254          DO  j = 1, radius_classes
     255             DO  i =  1, radius_classes
     256                ckernel(pstart+i-1,pstart+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j)
    130257             ENDDO
    131 
    132           ELSE
    133 
    134              CALL turbsd
    135              CALL turb_enhance_eff
    136              CALL effic
    137 
    138              DO  j = pstart, pend, 1
    139                 DO  i =  pstart, pend, 1
    140                    kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
    141                 ENDDO
     258          ENDDO
     259
     260          DEALLOCATE( gck, ecf )
     261
     262       ELSE
     263!
     264!--       Call routines to calculate efficiencies for the Hall kernel
     265          CALL fallg
     266          CALL effic
     267
     268          DO  j = 1, radius_classes
     269             DO  i =  1, radius_classes
     270                ckernel(pstart+i-1,pstart+j-1,1) = pi *                       &
     271                                          ( radclass(j) + radclass(i) )**2    &
     272                                          * ec(i,j) * ABS( winf(j) - winf(i) )
    142273             ENDDO
    143 
    144           ENDIF
    145 
    146           DEALLOCATE(gck, ecf)
    147 
    148        ELSE
    149 
    150 !          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
    151           CALL fallg
    152 !          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
    153 !          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
    154           CALL effic
    155 !          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
    156 
    157           DO  j = pstart, pend
    158              DO  i =  pstart, pend
    159                 kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
    160                                      particles(i)%radius * 100.0 )**2 &
    161                                  * ec(i,j) * ABS( winf(j) - winf(i) )
    162              ENDDO
    163274          ENDDO
    164275
    165276       ENDIF
    166277
    167        DEALLOCATE( ec, winf )
    168 
    169 !       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
    170 
    171     END SUBROUTINE colker
    172 
    173 !------------------------------------------------------------------------------!
    174 ! SUBROUTINE for calculation of w, g and gck
     278       ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
     279
     280       DEALLOCATE( ec, radclass, winf )
     281
     282    END SUBROUTINE recalculate_kernel
     283
     284
     285!------------------------------------------------------------------------------!
     286! Calculation of gck
     287! This is from Aayala 2008b, page 37ff.
     288! Necessary input parameters: water density, radii of droplets, air density,
     289! air viscosity, turbulent dissipation rate, taylor microscale reynolds number,
     290! gravitational acceleration  --> to be replaced by PALM parameters
    175291!------------------------------------------------------------------------------!
    176292    SUBROUTINE turbsd
    177 ! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
    178 ! of droplets, air density, air viscosity, turbulent dissipation rate,
    179 ! taylor microscale reynolds number, gravitational acceleration
    180293
    181294       USE constants
     
    186299       IMPLICIT NONE
    187300
    188        REAL :: Relamda, &
    189                Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
    190                bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
    191                vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
    192                wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
    193 
    194        REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
    195 
    196        REAL, DIMENSION(pstart:pend) :: St, tau
     301       INTEGER ::  i, j
    197302
    198303       LOGICAL, SAVE ::  first = .TRUE.
    199304
    200        INTEGER :: i, j
    201 
    202 !
    203 !--   Initial assignment of constants
     305       REAL ::  ao, ao_gr, bbb, be, b1, b2, ccc, c1, c1_gr, c2, d1, d2, eta, &
     306                e1, e2, fao_gr, fr, grfin, lambda, lambda_re, lf, rc, rrp,   &
     307                sst, tauk, tl, t2, tt, t1, vk, vrms1xy, vrms2xy, v1, v1v2xy, &
     308                v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z
     309
     310       REAL, SAVE ::  airdens, airvisc, anu, gravity, waterdens
     311
     312       REAL, DIMENSION(1:radius_classes) ::  st, tau
     313
     314
     315!
     316!--    Initial assignment of constants
    204317       IF ( first )  THEN
    205318
    206319          first = .FALSE.
    207           airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
    208           airdens   = 1.2250   !air density in mg/cm**3
    209           waterdens = 1000.0   !water density in mg/cm**3
    210           gravity   = 980.6650 !in cm/s**2
     320          airvisc   = 0.1818     ! dynamic viscosity in mg/cm*s
     321          airdens   = 1.2250     ! air density in mg/cm**3
     322          waterdens = 1000.0     ! water density in mg/cm**3
     323          gravity   = 980.6650   ! in cm/s**2
    211324          anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
    212325
    213326       ENDIF   
    214327
    215        Relamda = urms**2*sqrt(15.0/epsilon/anu)
    216 
    217        Tl = urms**2/epsilon       !in s
    218        Lf = 0.5 * (urms**3)/epsilon !in cm
    219 
    220        tauk = (anu/epsilon)**0.5       !in s
    221        eta  = (anu**3/epsilon)**0.25  !in cm
    222        vk   = eta/tauk                   !in cm/s
    223 
    224        ao = (11.+7.*Relamda)/(205.+Relamda)
    225 
    226        lambda = urms * sqrt(15.*anu/epsilon)     !in cm
    227 
    228        tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
    229 
    230        CALL fallg !gives winf in cm/s
    231 
    232        DO i = pstart, pend
    233           tau(i) = winf(i)/gravity    !in s
    234           St(i)  = tau(i)/tauk
     328       lambda    = urms * SQRT( 15.0 * anu / epsilon )     ! in cm
     329       lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu )
     330       tl        = urms**2 / epsilon                       ! in s
     331       lf        = 0.5 * urms**3 / epsilon                 ! in cm
     332       tauk      = SQRT( anu / epsilon )                   ! in s
     333       eta       = ( anu**3 / epsilon )**0.25              ! in cm
     334       vk        = eta / tauk                              ! in cm/s
     335
     336       ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
     337       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
     338
     339       CALL fallg    ! gives winf in cm/s
     340
     341       DO  i = 1, radius_classes
     342          tau(i) = winf(i) / gravity    ! in s
     343          st(i)  = tau(i) / tauk
    235344       ENDDO
    236345
    237 !***** TO CALCULATE wr ********************************
    238 !from Aayala 2008b, page 38f
    239 
    240        z  = tt/Tl
    241        be = sqrt(2.0)*lambda/Lf
    242 
    243        bbb = sqrt(1.0-2.0*be**2)
    244        d1 = (1.+bbb)/2.0/bbb
    245        e1 = Lf*(1.0+bbb)/2.0   !in cm
    246        d2 = (1.0-bbb)/2.0/bbb
    247        e2 = Lf*(1.0-bbb)/2.0   !in cm
    248 
    249        ccc = sqrt(1.0-2.0*z**2)
    250        b1 = (1.+ccc)/2./ccc
    251        c1 = Tl*(1.+ccc)/2.   !in s
    252        b2 = (1.-ccc)/2./ccc
    253        c2 = Tl*(1.-ccc)/2.   !in s
    254 
    255        DO i = pstart, pend
    256           v1 = winf(i)         !in cm/s
    257           t1 = tau(i)         !in s
    258 
    259           DO j = pstart,i
    260              rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
    261              v2  = winf(j)         !in cm/s
    262              t2  = tau(j)         !in s
    263 
    264              v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
    265                      - b1*d2*PHI(c1,e2,v1,t1) &
    266                      - b2*d1*PHI(c2,e1,v1,t1) &
    267                      + b2*d2*PHI(c2,e2,v1,t1)
    268              v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
    269 
    270              vrms1xy= sqrt(v1xysq)             !in cm/s
    271 
    272              v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
    273                      - b1*d2*PHI(c1,e2,v2,t2) &
    274                      - b2*d1*PHI(c2,e1,v2,t2) &
    275                      + b2*d2*PHI(c2,e2,v2,t2)
    276              v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
    277 
    278               vrms2xy= sqrt(v2xysq)             !in cm/s
    279 
    280              IF(winf(i).ge.winf(j)) THEN
     346!
     347!--    Calculate wr (from Aayala 2008b, page 38f)
     348       z   = tt / tl
     349       be  = SQRT( 2.0 ) * lambda / lf
     350       bbb = SQRT( 1.0 - 2.0 * be**2 )
     351       d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
     352       e1  = lf * ( 1.0 + bbb ) * 0.5   ! in cm
     353       d2  = ( 1.0 - bbb ) * 0.5 / bbb
     354       e2  = lf * ( 1.0 - bbb ) * 0.5   ! in cm
     355       ccc = SQRT( 1.0 - 2.0 * z**2 )
     356       b1  = ( 1.0 + ccc ) * 0.5 / ccc
     357       c1  = tl * ( 1.0 + ccc ) * 0.5   ! in s
     358       b2  = ( 1.0 - ccc ) * 0.5 / ccc
     359       c2  = tl * ( 1.0 - ccc ) * 0.5   ! in s
     360
     361       DO  i = 1, radius_classes
     362
     363          v1 = winf(i)        ! in cm/s
     364          t1 = tau(i)         ! in s
     365
     366          DO  j = 1, i
     367             rrp = radclass(i) + radclass(j)               ! radius in cm
     368             v2  = winf(j)                                 ! in cm/s
     369             t2  = tau(j)                                  ! in s
     370
     371             v1xysq  = b1 * d1 * phi(c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) &
     372                     - b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1)
     373             v1xysq  = v1xysq * urms**2 / t1                ! in cm**2/s**2
     374             vrms1xy = SQRT( v1xysq )                       ! in cm/s
     375
     376             v2xysq  = b1 * d1 * phi(c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) &
     377                     - b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2)
     378             v2xysq  = v2xysq * urms**2 / t2                ! in cm**2/s**2
     379             vrms2xy = SQRT( v2xysq )                       ! in cm/s
     380
     381             IF ( winf(i) >= winf(j) )  THEN
    281382                v1 = winf(i)
    282383                t1 = tau(i)
     
    290391             ENDIF
    291392
    292              v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
    293                      - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
    294                      - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
    295                      + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
    296              fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
    297              v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
    298 
    299              wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
    300 
    301               IF (wrtur2xy.le.0.0) wrtur2xy=0.0
    302 
    303               wrgrav2=pi/8.*(winf(j)-winf(i))**2
    304 
    305               wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
    306 
    307 
    308 !***** TO CALCULATE gr ********************************
    309 
    310              IF(St(j).gt.St(i)) THEN
    311                 SSt = St(j)
     393             v1v2xy   =  b1 * d1 * zhi(c1,e1,v1,t1,v2,t2) - &
     394                         b1 * d2 * zhi(c1,e2,v1,t1,v2,t2) - &
     395                         b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + &
     396                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
     397             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
     398             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in cm**2/s**2
     399             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in cm**2/s**2
     400             IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
     401             wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
     402             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in cm/s
     403
     404!
     405!--          Calculate gr
     406             IF ( st(j) > st(i) )  THEN
     407                sst = st(j)
    312408             ELSE
    313                 SSt = St(i)
     409                sst = st(i)
    314410             ENDIF
    315411
    316              XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
    317                   -4.2942*SSt**2 + 5.3406*SSt
    318 
    319              IF(XX.le.0.0) XX = 0.0
    320 
    321              YY = 0.1886*exp(20.306/Relamda)
    322 
    323              c1_gr = XX/(gravity/(vk/tauk))**YY
    324 
    325              ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
    326              fao_gr = 20.115 * (ao_gr/Relamda)**0.5
    327              rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
    328 
    329              grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
    330              IF (grFIN.lt.1.0) grFIN = 1.0
    331 
    332              gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
     412             xx = -0.1988 * sst**4 + 1.5275 * sst**3 - 4.2942 * sst**2 + &
     413                   5.3406 * sst
     414             IF ( xx < 0.0 )  xx = 0.0
     415             yy = 0.1886 * EXP( 20.306 / lambda_re )
     416
     417             c1_gr  =  xx / ( gravity / vk * tauk )**yy
     418
     419             ao_gr  = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2
     420             fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
     421             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
     422
     423             grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5 )
     424             IF ( grfin < 1.0 )  grfin = 1.0
     425
     426             gck(i,j) = 2.0 * pi * rrp**2 * wrfin * grfin           ! in cm**3/s
    333427             gck(j,i) = gck(i,j)
    334428
     
    336430       ENDDO
    337431
    338     END SUBROUTINE TurbSD
    339 
    340 !------------------------------------------------------------------------------!
    341 ! PHI as a function
    342 !------------------------------------------------------------------------------!
    343   REAL FUNCTION PHI(a,b,vsett,tau0)
     432    END SUBROUTINE turbsd
     433
     434
     435!------------------------------------------------------------------------------!
     436! phi as a function
     437!------------------------------------------------------------------------------!
     438    REAL FUNCTION phi( a, b, vsett, tau0 )
    344439
    345440       IMPLICIT NONE
    346441
    347        REAL :: a, aa1, b, vsett, tau0
    348 
    349        aa1 = 1./tau0 + 1./a + vsett/b
    350 
    351        PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
    352 
    353     END FUNCTION PHI
    354 
    355 !------------------------------------------------------------------------------!
    356 ! ZETA as a function
    357 !------------------------------------------------------------------------------!
    358   REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
     442       REAL ::  a, aa1, b, tau0, vsett
     443
     444       aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
     445       phi = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
     446
     447    END FUNCTION phi
     448
     449
     450!------------------------------------------------------------------------------!
     451! zeta as a function
     452!------------------------------------------------------------------------------!
     453    REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
    359454
    360455       IMPLICIT NONE
    361456
    362        REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
    363 
    364         aa1 = vsett2/b - 1./tau2 - 1./a
    365         aa2 = vsett1/b + 1./tau1 + 1./a
    366         aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
    367         aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
    368         aa5 = vsett2/b + 1./tau2 + 1./a
    369         aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
    370         ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
    371              + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
    372              + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
    373                                                            * 1./2./b/aa3             ! in s**2
    374 
    375     END FUNCTION ZHI
    376 
    377 !------------------------------------------------------------------------------!
    378 ! SUBROUTINE for calculation of terminal velocity winf
    379 !------------------------------------------------------------------------------!
    380    SUBROUTINE fallg
     457       REAL ::  a, aa1, aa2, aa3, aa4, aa5, aa6, b, tau1, tau2, vsett1, vsett2
     458
     459       aa1 = vsett2 / b - 1.0 / tau2 - 1.0 / a
     460       aa2 = vsett1 / b + 1.0 / tau1 + 1.0 / a
     461       aa3 = ( vsett1 - vsett2 ) / b + 1.0 / tau1 + 1.0 / tau2
     462       aa4 = ( vsett2 / b )**2 - ( 1.0 / tau2 + 1.0 / a )**2
     463       aa5 = vsett2 / b + 1.0 / tau2 + 1.0 / a
     464       aa6 = 1.0 / tau1 - 1.0 / a + ( 1.0 / tau2 + 1.0 / a) * vsett1 / vsett2
     465       zhi = (1.0 / aa1 - 1.0 / aa2 ) * ( vsett1 - vsett2 ) * 0.5 / b / aa3**2 &
     466           + (4.0 / aa4 - 1.0 / aa5**2 - 1.0 / aa1**2 ) * vsett2 * 0.5 / b /aa6&
     467           + (2.0 * ( b / aa2 - b / aa1 ) - vsett1 / aa2**2 + vsett2 / aa1**2 )&
     468           * 0.5 / b / aa3      ! in s**2
     469
     470    END FUNCTION zhi
     471
     472
     473!------------------------------------------------------------------------------!
     474! Calculation of terminal velocity winf
     475!------------------------------------------------------------------------------!
     476    SUBROUTINE fallg
    381477
    382478       USE constants
     
    385481       USE arrays_3d
    386482
    387       IMPLICIT NONE
    388 
    389       INTEGER :: i, j
    390 
    391       LOGICAL, SAVE ::  first = .TRUE.
    392 
    393       REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
    394 
    395       REAL :: bond, x, xrey, y
    396 
    397       REAL, DIMENSION(1:7), SAVE  :: b
    398       REAL, DIMENSION(1:6), SAVE  :: c
    399 
    400 !
    401 !--   Initial assignment of constants
    402       IF ( first )  THEN
     483       IMPLICIT NONE
     484
     485       INTEGER ::  i, j
     486
     487       LOGICAL, SAVE ::  first = .TRUE.
     488
     489       REAL, SAVE ::  cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &
     490                      t0, xlamb
     491
     492       REAL ::  bond, x, xrey, y
     493
     494       REAL, DIMENSION(1:7), SAVE  ::  b
     495       REAL, DIMENSION(1:6), SAVE  ::  c
     496
     497!
     498!--    Initial assignment of constants
     499       IF ( first )  THEN
     500
     501          first = .FALSE.
     502          b = (/  -0.318657E1,  0.992696E0, -0.153193E-2, -0.987059E-3, &
     503                 -0.578878E-3, 0.855176E-4, -0.327815E-5 /)
     504          c = (/  -0.500015E1,  0.523778E1,  -0.204914E1,   0.475294E0, &
     505                 -0.542819E-1, 0.238449E-2 /)
     506
     507          eta   = 1.818E-4          ! in poise = g/(cm s)
     508          xlamb = 6.62E-6           ! in cm
     509          rhow  = 1.0               ! in g/cm**3
     510          rhoa  = 1.225E-3          ! in g/cm**3
     511          grav  = 980.665           ! in cm/s**2
     512          cunh  = 1.257 * xlamb     ! in cm
     513          t0    = 273.15            ! in K
     514          sigma = 76.1 - 0.155 * ( 293.15 - t0 )              ! in N/m = g/s**2
     515          stok  = 2.0  * grav * ( rhow - rhoa ) / ( 9.0 * eta ) ! in 1/(cm s)
     516          stb   = 32.0 * rhoa * ( rhow - rhoa) * grav / (3.0 * eta * eta)
     517                                                                ! in 1/cm**3
     518          phy   = sigma**3 * rhoa**2 / ( eta**4 * grav * ( rhow - rhoa ) )
     519          py    = phy**( 1.0 / 6.0 )
     520
     521       ENDIF
     522
     523       DO  j = 1, radius_classes
     524
     525          IF ( radclass(j) <= 1.0E-3 ) THEN
     526
     527             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s
     528
     529          ELSEIF ( radclass(j) > 1.0E-3  .AND.  radclass(j) <= 5.35E-2 )  THEN
     530
     531             x = LOG( stb * radclass(j)**3 )
     532             y = 0.0
     533
     534             DO  i = 1, 7
     535                y = y + b(i) * x**(i-1)
     536             ENDDO
     537
     538             xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
     539             winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )      ! in cm/s
     540
     541          ELSEIF ( radclass(j) > 5.35E-2 )  THEN
     542
     543             IF ( radclass(j) > 0.35 )  THEN
     544                bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma
     545             ELSE
     546                bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma
     547             ENDIF
     548
     549             x = LOG( 16.0 * bond * py / 3.0 )
     550             y = 0.0
     551
     552             DO  i = 1, 6
     553                y = y + c(i) * x**(i-1)
     554             ENDDO
     555
     556             xrey = py * EXP( y )
     557
     558             IF ( radclass(j) > 0.35 )  THEN
     559                winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 )           ! in cm/s
     560             ELSE
     561                winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )    ! in cm/s
     562             ENDIF
     563
     564          ENDIF
     565
     566       ENDDO
     567
     568    END SUBROUTINE fallg
     569
     570
     571!------------------------------------------------------------------------------!
     572! Calculation of collision efficencies for the Hall kernel
     573!------------------------------------------------------------------------------!
     574    SUBROUTINE effic
     575
     576       USE arrays_3d
     577       USE cloud_parameters
     578       USE constants
     579       USE particle_attributes
     580
     581       IMPLICIT NONE
     582
     583       INTEGER ::  i, iq, ir, j, k, kk
     584
     585       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
     586
     587       LOGICAL, SAVE ::  first = .TRUE.
     588
     589       REAL ::  ek, particle_radius, pp, qq, rq
     590
     591       REAL, DIMENSION(1:21), SAVE ::  rat
     592       REAL, DIMENSION(1:15), SAVE ::  r0
     593       REAL, DIMENSION(1:15,1:21), SAVE ::  ecoll
     594
     595!
     596!--    Initial assignment of constants
     597       IF ( first )  THEN
    403598
    404599         first = .FALSE.
    405          b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
    406              -0.578878e-3,0.855176e-4,-0.327815e-5/)
    407          c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
    408              -0.542819e-1, 0.238449e-2/)
    409 
    410          eta   = 1.818e-4     !in poise = g/(cm s)
    411          xlamb = 6.62e-6    !in cm
    412 
    413          rhow = 1.0       !in g/cm**3
    414          rhoa = 1.225e-3    !in g/cm**3
    415 
    416          grav = 980.665     !in cm/s**2
    417          cunh = 1.257 * xlamb !in cm
    418          t0   = 273.15        !in K
    419          sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
    420          stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
    421          stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
    422          phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
    423          py   = phy**(1.0/6.0)
    424 
    425       ENDIF
    426 
    427 !particle radius has to be in cm
    428       DO j = pstart, pend
    429 
    430          IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
    431 
    432             winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
    433 
    434          ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
    435 
    436             x = log(stb*(particles(j)%radius*100.0)**3)
    437             y = 0.0
    438 
    439             DO i = 1, 7
    440                y = y + b(i) * (x**(i-1))
    441             ENDDO
    442 
    443             xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
    444             winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
    445 
    446          ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
    447 
    448             IF (particles(j)%radius*100.0.gt.0.35)  THEN
    449                bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
    450             ELSE
    451                bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
    452             ENDIF
    453 
    454             x = log(16.0*bond*py/3.0)
    455             y = 0.0
    456 
    457             DO i = 1, 6
    458                y = y + c(i) * (x**(i-1))
    459             ENDDO
    460 
    461             xrey = py*exp(y)
    462 
    463             IF (particles(j)%radius*100.0 .gt.0.35)  THEN
    464               winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
    465             ELSE
    466                winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
    467             ENDIF
    468 
    469          ENDIF
    470       ENDDO
    471       RETURN
    472    END SUBROUTINE fallg
    473 
    474 !------------------------------------------------------------------------------!
    475 ! SUBROUTINE for calculation of collision efficencies
    476 !------------------------------------------------------------------------------!
    477 
    478    SUBROUTINE effic
    479 
    480       USE arrays_3d
    481       USE constants
    482       USE cloud_parameters
    483       USE particle_attributes
    484 
    485 !collision efficiencies of hall kernel
    486       IMPLICIT NONE
    487 
    488       INTEGER :: i, ir, iq, j, k, kk
    489 
    490       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
    491 
    492       LOGICAL, SAVE ::  first = .TRUE.
    493 
    494       REAL ::  ek, particle_radius, pp, qq, rq
    495 
    496       REAL, DIMENSION(1:21), SAVE :: rat
    497       REAL, DIMENSION(1:15), SAVE :: r0
    498       REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
    499 
    500 !
    501 !--   Initial assignment of constants
    502       IF ( first )  THEN
    503 
    504          first = .FALSE.
    505          r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
    506                  300./)
    507          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,  &
    508                  0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
    509 
    510          ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
    511                         0.001,0.001,0.001,0.001,0.001,0.001/)
    512          ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
    513                         0.050,0.200,0.500,0.770,0.870,0.970/)
    514          ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
    515                         0.430,0.580,0.790,0.930,0.960,1.000/)
    516          ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
    517                         0.640,0.750,0.910,0.970,0.980,1.000/)
    518          ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
    519                         0.770,0.840,0.950,0.970,1.000,1.000/)
    520          ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
    521                         0.840,0.880,0.950,1.000,1.000,1.000/)
    522          ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
    523                         0.870,0.900,0.950,1.000,1.000,1.000/)
    524          ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
    525                         0.890,0.920,1.000,1.000,1.000,1.000/)
    526          ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
    527                         0.900,0.940,1.000,1.000,1.000,1.000/)
    528          ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
    529                         0.910,0.950,1.000,1.000,1.000,1.000/)
    530          ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
    531                         0.910,0.950,1.000,1.000,1.000,1.000/)
    532          ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
    533                         0.910,0.950,1.000,1.000,1.000,1.000/)
    534          ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
    535                         0.910,0.950,1.000,1.000,1.000,1.000/)
    536          ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
    537                         0.910,0.950,1.000,1.000,1.000,1.000/)
    538          ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
    539                         0.920,0.950,1.000,1.000,1.000,1.000/)
    540          ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
    541                         0.930,0.970,1.000,1.000,1.000,1.000/)
    542          ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
    543                         0.950,1.000,1.000,1.000,1.000,1.000/)
    544          ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
    545                         1.000,1.020,1.020,1.020,1.020,1.020/)
    546          ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
    547                         1.030,1.040,1.040,1.040,1.040,1.040/)
    548          ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
    549                         1.700,2.300,2.300,2.300,2.300,2.300/)
    550          ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
    551                         3.000,4.000,4.000,4.000,4.000,4.000/)
    552       ENDIF
    553 
    554 !
    555 !--   Calculate the radius class index of particles with respect to array r
    556       ALLOCATE( ira(pstart:pend) )
    557       DO  j = pstart, pend
    558          particle_radius = particles(j)%radius * 1.0E6
    559          DO  k = 1, 15
    560             IF ( particle_radius < r0(k) )  THEN
    561                ira(j) = k
    562                EXIT
    563             ENDIF
    564          ENDDO
    565          IF ( particle_radius >= r0(15) )  ira(j) = 16
    566       ENDDO
    567 
    568 !
    569 !--   Two-dimensional linear interpolation of the collision efficiency.
    570 !--   Radius has to be in µm
    571       DO  j = pstart, pend
    572          DO  i = pstart, j
    573 
    574             ir = ira(j)
    575 
    576             rq = particles(i)%radius / particles(j)%radius
    577 
    578 !            DO kk = 2, 21
    579 !               IF ( rq <= rat(kk) )  THEN
    580 !                  iq = kk
    581 !                  EXIT
    582 !               ENDIF
    583 !            ENDDO
    584 
    585             iq = INT( rq * 20 ) + 1
    586             iq = MAX(iq , 2)
    587 
    588             IF ( ir < 16 )  THEN
    589 
    590                IF ( ir >= 2 )  THEN
    591                   pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
    592                        ( r0(ir) - r0(ir-1) )
    593                   qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    594                   ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
    595                             + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
    596                             + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
    597                             + pp * qq * ecoll(ir,iq)
    598                ELSE
    599                   qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    600                   ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
    601                ENDIF
    602 
    603             ELSE
    604                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    605                ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
    606                ec(j,i) = MIN( ek, 1.0 )
    607             ENDIF
    608 
    609             ec(i,j) = ec(j,i)
    610             IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
    611 
    612          ENDDO
    613       ENDDO
    614 
    615       DEALLOCATE( ira )
    616 
    617    END SUBROUTINE effic
    618 
    619 !------------------------------------------------------------------------------!
    620 ! SUBROUTINE for calculation of enhancement factor collision efficencies
    621 !------------------------------------------------------------------------------!
    622    SUBROUTINE turb_enhance_eff
     600         r0  = (/ 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60., &
     601                  70.0, 100.0, 150.0, 200.0, 300.0 /)
     602         rat = (/ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, &
     603                  0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, &
     604                  1.00 /)
     605
     606         ecoll(:,1) = (/0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, &
     607                        0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001/)
     608         ecoll(:,2) = (/0.003, 0.003, 0.003, 0.004, 0.005, 0.005, 0.005, &
     609                        0.010, 0.100, 0.050, 0.200, 0.500, 0.770, 0.870, 0.970/)
     610         ecoll(:,3) = (/0.007, 0.007, 0.007, 0.008, 0.009, 0.010, 0.010, &
     611                        0.070, 0.400, 0.430, 0.580, 0.790, 0.930, 0.960, 1.000/)
     612         ecoll(:,4) = (/0.009, 0.009, 0.009, 0.012, 0.015, 0.010, 0.020, &
     613                        0.280, 0.600, 0.640, 0.750, 0.910, 0.970, 0.980, 1.000/)
     614         ecoll(:,5) = (/0.014, 0.014, 0.014, 0.015, 0.016, 0.030, 0.060, &
     615                        0.500, 0.700, 0.770, 0.840, 0.950, 0.970, 1.000, 1.000/)
     616         ecoll(:,6) = (/0.017, 0.017, 0.017, 0.020, 0.022, 0.060, 0.100, &
     617                        0.620, 0.780, 0.840, 0.880, 0.950, 1.000, 1.000, 1.000/)
     618         ecoll(:,7) = (/0.030, 0.030, 0.024, 0.022, 0.032, 0.062, 0.200, &
     619                        0.680, 0.830, 0.870, 0.900, 0.950, 1.000, 1.000, 1.000/)
     620         ecoll(:,8) = (/0.025, 0.025, 0.025, 0.036, 0.043, 0.130, 0.270, &
     621                        0.740, 0.860, 0.890, 0.920, 1.000, 1.000, 1.000, 1.000/)
     622         ecoll(:,9) = (/0.027, 0.027, 0.027, 0.040, 0.052, 0.200, 0.400, &
     623                        0.780, 0.880, 0.900, 0.940, 1.000, 1.000, 1.000, 1.000/)
     624         ecoll(:,10)= (/0.030, 0.030, 0.030, 0.047, 0.064, 0.250, 0.500, &
     625                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     626         ecoll(:,11)= (/0.040, 0.040, 0.033, 0.037, 0.068, 0.240, 0.550, &
     627                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     628         ecoll(:,12)= (/0.035, 0.035, 0.035, 0.055, 0.079, 0.290, 0.580, &
     629                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     630         ecoll(:,13)= (/0.037, 0.037, 0.037, 0.062, 0.082, 0.290, 0.590, &
     631                        0.780, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     632         ecoll(:,14)= (/0.037, 0.037, 0.037, 0.060, 0.080, 0.290, 0.580, &
     633                        0.770, 0.890, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     634         ecoll(:,15)= (/0.037, 0.037, 0.037, 0.041, 0.075, 0.250, 0.540, &
     635                        0.760, 0.880, 0.920, 0.950, 1.000, 1.000, 1.000, 1.000/)
     636         ecoll(:,16)= (/0.037, 0.037, 0.037, 0.052, 0.067, 0.250, 0.510, &
     637                        0.770, 0.880, 0.930, 0.970, 1.000, 1.000, 1.000, 1.000/)
     638         ecoll(:,17)= (/0.037, 0.037, 0.037, 0.047, 0.057, 0.250, 0.490, &
     639                        0.770, 0.890, 0.950, 1.000, 1.000, 1.000, 1.000, 1.000/)
     640         ecoll(:,18)= (/0.036, 0.036, 0.036, 0.042, 0.048, 0.230, 0.470, &
     641                        0.780, 0.920, 1.000, 1.020, 1.020, 1.020, 1.020, 1.020/)
     642         ecoll(:,19)= (/0.040, 0.040, 0.035, 0.033, 0.040, 0.112, 0.450, &
     643                        0.790, 1.010, 1.030, 1.040, 1.040, 1.040, 1.040, 1.040/)
     644         ecoll(:,20)= (/0.033, 0.033, 0.033, 0.033, 0.033, 0.119, 0.470, &
     645                        0.950, 1.300, 1.700, 2.300, 2.300, 2.300, 2.300, 2.300/)
     646         ecoll(:,21)= (/0.027, 0.027, 0.027, 0.027, 0.027, 0.125, 0.520, &
     647                        1.400, 2.300, 3.000, 4.000, 4.000, 4.000, 4.000, 4.000/)
     648       ENDIF
     649
     650!
     651!--    Calculate the radius class index of particles with respect to array r
     652       ALLOCATE( ira(1:radius_classes) )
     653       DO  j = 1, radius_classes
     654          particle_radius = radclass(j) * 1.0E4   ! in microm
     655          DO  k = 1, 15
     656             IF ( particle_radius < r0(k) )  THEN
     657                ira(j) = k
     658                EXIT
     659             ENDIF
     660          ENDDO
     661          IF ( particle_radius >= r0(15) )  ira(j) = 16
     662       ENDDO
     663
     664!
     665!--    Two-dimensional linear interpolation of the collision efficiency.
     666!--    Radius has to be in µm
     667       DO  j = 1, radius_classes
     668          DO  i = 1, j
     669
     670             ir = ira(j)
     671             rq = radclass(i) / radclass(j)
     672             iq = INT( rq * 20 ) + 1
     673             iq = MAX( iq , 2)
     674
     675             IF ( ir < 16 )  THEN
     676                IF ( ir >= 2 )  THEN
     677                   pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / &
     678                        ( r0(ir) - r0(ir-1) )
     679                   qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     680                   ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)  &
     681                             + pp * ( 1.0-qq ) * ecoll(ir,iq-1)          &
     682                             + qq * ( 1.0-pp ) * ecoll(ir-1,iq)          &
     683                             + pp * qq * ecoll(ir,iq)
     684                ELSE
     685                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     686                   ec(j,i) = (1.0-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
     687                ENDIF
     688             ELSE
     689                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     690                ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
     691                ec(j,i) = MIN( ek, 1.0 )
     692             ENDIF
     693
     694             ec(i,j) = ec(j,i)
     695             IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
     696
     697          ENDDO
     698       ENDDO
     699
     700       DEALLOCATE( ira )
     701
     702    END SUBROUTINE effic
     703
     704
     705!------------------------------------------------------------------------------!
     706! Calculation of enhancement factor for collision efficencies due to turbulence
     707!------------------------------------------------------------------------------!
     708    SUBROUTINE turb_enhance_eff
    623709
    624710       USE constants
     
    627713       USE arrays_3d
    628714
    629       IMPLICIT NONE
    630 
    631       INTEGER :: i, ik, ir, iq, j, k, kk
    632 
    633       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
    634 
    635       REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
    636 
    637       LOGICAL, SAVE ::  first = .TRUE.
    638 
    639       REAL, DIMENSION(1:11), SAVE :: rat
    640       REAL, DIMENSION(1:7), SAVE  :: r0
    641       REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
    642 
    643 !
    644 !--   Initial assignment of constants
    645       IF ( first )  THEN
    646 
    647          first = .FALSE.
    648 
    649          r0  = (/10., 20., 30.,40., 50., 60.,100./)
    650          rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
    651 
    652 ! 100 cm^2/s^3
    653          ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
    654          ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
    655          ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
    656          ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
    657          ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
    658          ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
    659          ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
    660          ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
    661          ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
    662          ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
    663          ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
    664 
    665 ! 400 cm^2/s^3
    666          ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
    667          ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
    668          ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
    669          ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
    670          ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
    671          ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
    672          ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
    673          ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
    674          ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
    675          ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
    676          ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
    677 
    678       ENDIF
    679 
    680 !
    681 !--   Calculate the radius class index of particles with respect to array r
    682       ALLOCATE( ira(pstart:pend) )
    683 
    684       DO  j = pstart, pend
    685          particle_radius = particles(j)%radius * 1.0E6
    686          DO  k = 1, 7
    687             IF ( particle_radius < r0(k) )  THEN
    688                ira(j) = k
    689                EXIT
    690             ENDIF
    691          ENDDO
    692          IF ( particle_radius >= r0(7) )  ira(j) = 8
    693       ENDDO
    694 
    695 ! two-dimensional linear interpolation of the collision efficiency
    696       DO j =  pstart, pend
    697          DO i = pstart, j
    698 
    699             ir = ira(j)
    700 
    701             rq = particles(i)%radius/particles(j)%radius
    702 
    703             DO kk = 2, 11
    704                IF ( rq <= rat(kk) )  THEN
    705                   iq = kk
    706                   EXIT
    707                ENDIF
    708             ENDDO
    709 
    710 ! 0  cm2/s3
    711             y1 = 1.0
    712 ! 100 cm2/s3, 400 cm2/s3
    713             IF (ir.lt.8) THEN
    714                IF (ir.ge.2) THEN
    715                   pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
    716                   qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    717                   y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
    718                             pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
    719                             qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
    720                             pp*qq*ecoll_100(ir,iq)
    721                   y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
    722                             pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
    723                             qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
    724                             pp*qq*ecoll_400(ir,iq)
    725                ELSE
    726                   qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    727                   y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
    728                   y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
    729                ENDIF
    730             ELSE
    731             qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    732             y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
    733             y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
    734          ENDIF
    735 ! linear interpolation
    736 ! dissipation rate in cm2/s3
    737          x1 = 0.0
    738          x2 = 100.0
    739          x3 = 400.0
    740 
    741          IF (epsilon.le.100.) THEN
    742             ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
    743                     +  (epsilon-0.)/(100.-0.) * y2
    744          ELSE IF(epsilon.le.600.)THEN
    745             ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
    746                     +  (epsilon-100.)/(400.-100.) * y3
    747 
    748          ELSE
    749             ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
    750                     +  (600.-100.)/(400.-100.) * y3
    751          ENDIF
    752 
    753          IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
    754 
    755          ecf(i,j)=ecf(j,i)
    756       ENDDO
    757    ENDDO
    758 
    759    RETURN
    760    END SUBROUTINE turb_enhance_eff
     715       IMPLICIT NONE
     716
     717       INTEGER :: i, ik, iq, ir, j, k, kk
     718
     719       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
     720
     721       REAL ::  particle_radius, pp, qq, rq, x1, x2, x3, y1, y2, y3
     722
     723       LOGICAL, SAVE ::  first = .TRUE.
     724
     725       REAL, DIMENSION(1:11), SAVE ::  rat
     726       REAL, DIMENSION(1:7), SAVE  ::  r0
     727       REAL, DIMENSION(1:7,1:11), SAVE ::  ecoll_100, ecoll_400
     728
     729!
     730!--    Initial assignment of constants
     731       IF ( first )  THEN
     732
     733          first = .FALSE.
     734
     735          r0  = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 100.0 /)
     736          rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
     737!
     738!--       In 100 cm^2/s^3
     739          ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
     740          ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
     741          ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
     742          ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
     743          ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
     744          ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
     745          ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
     746          ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
     747          ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
     748          ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
     749          ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
     750!
     751!--       In 400 cm^2/s^3
     752          ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
     753          ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
     754          ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
     755          ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
     756          ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
     757          ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
     758          ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
     759          ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
     760          ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
     761          ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
     762          ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
     763
     764       ENDIF
     765
     766!
     767!--    Calculate the radius class index of particles with respect to array r0
     768       ALLOCATE( ira(1:radius_classes) )
     769
     770       DO  j = 1, radius_classes
     771          particle_radius = radclass(j) * 1.0E4    ! in microm
     772          DO  k = 1, 7
     773             IF ( particle_radius < r0(k) )  THEN
     774                ira(j) = k
     775                EXIT
     776             ENDIF
     777          ENDDO
     778          IF ( particle_radius >= r0(7) )  ira(j) = 8
     779       ENDDO
     780
     781!
     782!--    Two-dimensional linear interpolation of the collision efficiencies
     783       DO  j =  1, radius_classes
     784          DO  i = 1, j
     785
     786             ir = ira(j)
     787             rq = radclass(i) / radclass(j)
     788
     789             DO  kk = 2, 11
     790                IF ( rq <= rat(kk) )  THEN
     791                   iq = kk
     792                   EXIT
     793                ENDIF
     794             ENDDO
     795
     796             y1 = 1.0      ! 0  cm2/s3
     797!
     798!--          100 cm2/s3, 400 cm2/s3
     799             IF ( ir < 8 )  THEN
     800                IF ( ir >= 2 )  THEN
     801                   pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
     802                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     803                   y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) +  &
     804                                pp * ( 1.0-qq ) * ecoll_100(ir,iq-1)   +  &
     805                                qq * ( 10.-pp ) * ecoll_100(ir-1,iq)   +  &
     806                                pp * qq         * ecoll_100(ir,iq)
     807                   y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) +  &
     808                                pp * ( 1.0-qq ) * ecoll_400(ir,iq-1)   +  &
     809                                qq * ( 1.0-pp ) * ecoll_400(ir-1,iq)   +  &
     810                                pp * qq         * ecoll_400(ir,iq)
     811                ELSE
     812                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     813                   y2 = ( 1.0-qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq)
     814                   y3 = ( 1.0-qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq)
     815                ENDIF
     816             ELSE
     817                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     818                y2 = ( 1.0-qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
     819                y3 = ( 1.0-qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
     820             ENDIF
     821!
     822!--          Linear interpolation of dissipation rate in cm2/s3
     823             IF ( epsilon <= 100.0 )  THEN
     824                ecf(j,i) = ( epsilon - 100.0 ) / (   0.0 - 100.0 ) * y1 &
     825                         + ( epsilon -   0.0 ) / ( 100.0 -   0.0 ) * y2
     826             ELSEIF ( epsilon <= 600.0 )  THEN
     827                ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
     828                         + ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3
     829             ELSE
     830                ecf(j,i) = (   600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
     831                         + (   600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3
     832             ENDIF
     833
     834             IF ( ecf(j,i) < 1.0 )  ecf(j,i) = 1.0
     835
     836             ecf(i,j) = ecf(j,i)
     837
     838          ENDDO
     839       ENDDO
     840
     841    END SUBROUTINE turb_enhance_eff
    761842
    762843 END MODULE lpm_collision_kernels_mod
  • palm/trunk/SOURCE/modules.f90

    r826 r828  
    44! Current revisions:
    55! -----------------
    6 !
     6! +dissipation_classes, radius_classes, use_kernel_tables,
     7! particle feature color renamed class
    78!
    89! Former revisions:
     
    11951196    INTEGER ::  mpi_particle_type
    11961197#endif
    1197     INTEGER ::  ibc_par_lr, ibc_par_ns, ibc_par_b, ibc_par_t,                  &
    1198                 iran_part = -1234567, maximum_number_of_particles = 1000,      &
     1198    INTEGER ::  dissipation_classes = 10, ibc_par_lr, ibc_par_ns, ibc_par_b,   &
     1199                ibc_par_t, iran_part = -1234567,                               &
     1200                maximum_number_of_particles = 1000,                            &
    11991201                maximum_number_of_tailpoints = 100,                            &
    12001202                maximum_number_of_tails = 0,                                   &
     
    12031205                number_of_initial_tails = 0, offset_ocean_nzt = 0,             &
    12041206                offset_ocean_nzt_m1 = 0, particles_per_point = 1,              &
    1205                 particle_file_count = 0, skip_particles_for_tail = 100,        &
    1206                 sort_count = 0, total_number_of_particles,                     &
    1207                 total_number_of_tails = 0
     1207                particle_file_count = 0, radius_classes = 20,                  &
     1208                skip_particles_for_tail = 100, sort_count = 0,                 &
     1209                total_number_of_particles, total_number_of_tails = 0
    12081210
    12091211    INTEGER, PARAMETER ::  max_number_of_particle_groups = 10
     
    12151217                particle_advection = .FALSE., random_start_position = .FALSE., &
    12161218                read_particles_from_restartfile = .TRUE.,                      &
    1217                 uniform_particles = .TRUE., use_particle_tails = .FALSE.,      &
    1218                 use_sgs_for_particles = .FALSE.,                              &
     1219                uniform_particles = .TRUE., use_kernel_tables = .FALSE.,       &
     1220                use_particle_tails = .FALSE., use_sgs_for_particles = .FALSE., &
    12191221                wang_kernel = .FALSE., write_particle_statistics = .FALSE.
    12201222
     
    12471249                   origin_z, radius, rvar1, rvar2, rvar3, speed_x, speed_y, &
    12481250                   speed_z, weight_factor, x, y, z
    1249        INTEGER ::  color, group, tailpoints, tail_id
     1251       INTEGER ::  class, group, tailpoints, tail_id
    12501252    END TYPE particle_type
    12511253
  • palm/trunk/SOURCE/package_parin.f90

    r826 r828  
    44! Current revisions:
    55! -----------------
    6 !
     6! +dissipation_classes, radius_classes in parpar
    77!
    88! Former revisions:
     
    7575
    7676    NAMELIST /particles_par/      bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,    &
    77                                   collision_kernel,                            &
    78                                   density_ratio, radius, dt_dopts,             &
     77                                  collision_kernel, density_ratio,             &
     78                                  dissipation_classes, dt_dopts,               &
    7979                                  dt_min_part, dt_prel, dt_sort_particles,     &
    8080                                  dt_write_particle_data,                      &
     
    8888                                  particle_advection_start,                    &
    8989                                  particle_maximum_age, pdx, pdy, pdz, psb,    &
    90                                   psl, psn, psr, pss, pst,                     &
    91                                   random_start_position,                       &
     90                                  psl, psn, psr, pss, pst, radius,             &
     91                                  radius_classes, random_start_position,       &
    9292                                  read_particles_from_restartfile,             &
    9393                                  skip_particles_for_tail,                     &
  • palm/trunk/SOURCE/set_particle_attributes.f90

    r623 r828  
    44! Current revisions:
    55! -----------------
    6 !
     6! particle feature color renamed class
    77!
    88! Former revisions:
     
    121121!
    122122!--       Number of available colors is defined in init_dvrp
    123           particles(n)%color = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
     123          particles(n)%class = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
    124124
    125125       ENDDO
     
    192192!
    193193!--       Number of available colors is defined in init_dvrp
    194           particles(n)%color = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
     194          particles(n)%class = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
    195195
    196196       ENDDO
     
    213213!
    214214!--       Number of available colors is defined in init_dvrp
    215           particles(n)%color = 1 + height * ( dvrp_colortable_entries_prt - 1 )
     215          particles(n)%class = 1 + height * ( dvrp_colortable_entries_prt - 1 )
    216216
    217217       ENDDO
Note: See TracChangeset for help on using the changeset viewer.