Changeset 1822

Timestamp:

Apr 7, 2016 7:49:42 AM (9 years ago)

Author:

hoffmann

Message:

changes in LPM and bulk cloud microphysics

Location:

palm/trunk/SOURCE

Files:

• Makefile (modified) (10 diffs)
• advec_ws.f90 (modified) (6 diffs)
• boundary_conds.f90 (modified) (7 diffs)
• calc_liquid_water_content.f90 (modified) (4 diffs)
• calc_precipitation.f90 (deleted)
• check_parameters.f90 (modified) (23 diffs)
• data_output_2d.f90 (modified) (3 diffs)
• data_output_3d.f90 (modified) (7 diffs)
• data_output_dvrp.f90 (modified) (4 diffs)
• flow_statistics.f90 (modified) (7 diffs)
• header.f90 (modified) (8 diffs)
• impact_of_latent_heat.f90 (deleted)
• init_3d_model.f90 (modified) (13 diffs)
• init_cloud_physics.f90 (modified) (3 diffs)
• init_masks.f90 (modified) (5 diffs)
• interaction_droplets_ptq.f90 (modified) (2 diffs)
• lpm.f90 (modified) (7 diffs)
• lpm_advec.f90 (modified) (10 diffs)
• lpm_boundary_conds.f90 (modified) (10 diffs)
• lpm_calc_liquid_water_content.f90 (modified) (2 diffs)
• lpm_collision_kernels.f90 (modified) (30 diffs)
• lpm_data_output_particles.f90 (modified) (4 diffs)
• lpm_droplet_collision.f90 (modified) (8 diffs)
• lpm_droplet_condensation.f90 (modified) (4 diffs)
• lpm_exchange_horiz.f90 (modified) (43 diffs)
• lpm_extend_tail_array.f90 (deleted)
• lpm_extend_tails.f90 (deleted)
• lpm_init.f90 (modified) (17 diffs)
• lpm_init_sgs_tke.f90 (modified) (2 diffs)
• lpm_pack_arrays.f90 (modified) (6 diffs)
• lpm_read_restart_file.f90 (modified) (7 diffs)
• lpm_release_set.f90 (deleted)
• lpm_set_attributes.f90 (modified) (4 diffs)
• lpm_write_exchange_statistics.f90 (modified) (2 diffs)
• lpm_write_restart_file.f90 (modified) (6 diffs)
• microphysics.f90 (modified) (45 diffs)
• mod_particle_attributes.f90 (modified) (5 diffs)
• modules.f90 (modified) (4 diffs)
• package_parin.f90 (modified) (5 diffs)
• prognostic_equations.f90 (modified) (22 diffs)
• surface_layer_fluxes.f90 (modified) (4 diffs)
• swap_timelevel.f90 (modified) (5 diffs)
• time_integration.f90 (modified) (5 diffs)
• user_actions.f90 (modified) (3 diffs)
• write_3d_binary.f90 (modified) (3 diffs)

palm/trunk/SOURCE/Makefile

@@ -20 +20 @@
 # Current revisions:
 # ------------------
-# 
+# Tails removed. lpm_release_set removed. calc_precipitation, impact_of_latent_heat
+# removed. 
 # 
 # Former revisions:
@@ -246 +247 @@
         advec_w_pw.f90 advec_w_up.f90 average_3d_data.f90 boundary_conds.f90 \
         buoyancy.f90 calc_liquid_water_content.f90 calc_mean_profile.f90 \
-        calc_precipitation.f90 calc_radiation.f90 \
+        calc_radiation.f90 \
         check_for_restart.f90 check_open.f90 check_parameters.f90 \
         close_file.f90 compute_vpt.f90 coriolis.f90 cpulog.f90 \
@@ -257 +258 @@
         disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \
         exchange_horiz_2d.f90 fft_xy.f90 flow_statistics.f90 \
-        global_min_max.f90 header.f90 impact_of_latent_heat.f90 \
+        global_min_max.f90 header.f90 \
         inflow_turbulence.f90 init_1d_model.f90 init_3d_model.f90 \
         init_advec.f90 init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 \
@@ -268 +269 @@
         lpm_data_output_particles.f90 lpm_droplet_collision.f90 \
         lpm_droplet_condensation.f90 lpm_exchange_horiz.f90 \
-        lpm_extend_tails.f90 \
-        lpm_extend_tail_array.f90 lpm_init.f90 lpm_init_sgs_tke.f90 \
-        lpm_pack_arrays.f90 lpm_read_restart_file.f90 lpm_release_set.f90 \
+        lpm_init.f90 lpm_init_sgs_tke.f90 \
+        lpm_pack_arrays.f90 lpm_read_restart_file.f90 \
         lpm_set_attributes.f90 \
         lpm_write_exchange_statistics.f90 lpm_write_restart_file.f90 \
@@ -349 +349 @@
 calc_mean_profile.o: modules.o mod_kinds.o
 calc_liquid_water_content.o: modules.o mod_kinds.o
-calc_precipitation.o: modules.o mod_kinds.o
 calc_radiation.o: modules.o mod_kinds.o
 check_for_restart.o: modules.o mod_kinds.o pmc_interface.o
@@ -395 +394 @@
    plant_canopy_model.o pmc_handle_communicator.o pmc_interface.o \
    radiation_model.o spectrum.o subsidence.o 
-impact_of_latent_heat.o: modules.o mod_kinds.o
 inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
 init_1d_model.o: modules.o mod_kinds.o
@@ -418 +416 @@
 local_tremain.o: modules.o cpulog.o mod_kinds.o
 local_tremain_ini.o: modules.o cpulog.o mod_kinds.o
-lpm.o: modules.o cpulog.o lpm_exchange_horiz.o lpm_pack_arrays.o mod_kinds.o \
-        mod_particle_attributes.o
+lpm.o: modules.o cpulog.o lpm_exchange_horiz.o lpm_init.o lpm_pack_arrays.o \
+        mod_kinds.o mod_particle_attributes.o
 lpm_advec.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
 lpm_boundary_conds.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
@@ -434 +432 @@
 lpm_exchange_horiz.o: modules.o cpulog.o lpm_pack_arrays.o mod_kinds.o \
         mod_particle_attributes.o netcdf_interface.o
-lpm_extend_tails.o: modules.o mod_kinds.o mod_particle_attributes.o
-lpm_extend_tail_array.o: modules.o mod_kinds.o mod_particle_attributes.o
 lpm_init.o: modules.o lpm_collision_kernels.o mod_kinds.o \
         netcdf_interface.o random_function.o mod_particle_attributes.o \
@@ -443 +439 @@
 lpm_read_restart_file.o: modules.o mod_kinds.o mod_particle_attributes.o \
         lpm_pack_arrays.o
-lpm_release_set.o: modules.o mod_kinds.o random_function.o \
-        mod_particle_attributes.o lpm_init.o  random_generator_parallel.o
 lpm_set_attributes.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
 lpm_write_exchange_statistics.o: modules.o mod_kinds.o mod_particle_attributes.o
@@ -480 +474 @@
 prognostic_equations.o: modules.o advec_s_pw.o advec_s_up.o advec_s_bc.o advec_u_pw.o \
         advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o \
-        advec_ws.o buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \
+        advec_ws.o buoyancy.o calc_radiation.o coriolis.o \
         cpulog.o diffusion_e.o diffusion_s.o diffusion_u.o diffusion_v.o diffusion_w.o \
-        eqn_state_seawater.o impact_of_latent_heat.o mod_kinds.o microphysics.o \
+        eqn_state_seawater.o mod_kinds.o microphysics.o \
         nudging.o plant_canopy_model.o production_e.o radiation_model.o \
         subsidence.o user_actions.o

palm/trunk/SOURCE/advec_ws.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! icloud_scheme removed, microphysics_seifert added
 ! 
 ! Former revisions:
@@ -222 +222 @@
 
        USE control_parameters,                                                 &
-           ONLY:  cloud_physics, humidity, icloud_scheme, loop_optimization,   &
-                  monotonic_adjustment, passive_scalar, precipitation, ocean,  &
-                  ws_scheme_mom, ws_scheme_sca
+           ONLY:  cloud_physics, humidity, loop_optimization,                  &
+                  monotonic_adjustment, passive_scalar, microphysics_seifert,  &
+                  ocean, ws_scheme_mom, ws_scheme_sca
 
        USE indices,                                                            &
@@ -274 +274 @@
           ENDIF
 
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
-               precipitation )  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
              ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
@@ -331 +330 @@
              ENDIF
 
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-                  precipitation )  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
                           diss_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
@@ -365 +363 @@
     
        USE control_parameters,                                                 &
-           ONLY:  cloud_physics, humidity, icloud_scheme, passive_scalar,      &
-                  precipitation, ocean, ws_scheme_mom, ws_scheme_sca
+           ONLY:  cloud_physics, humidity, passive_scalar, ocean,              &
+                  microphysics_seifert, ws_scheme_mom, ws_scheme_sca
 
        USE kinds
@@ -391 +389 @@
           sums_wspts_ws_l = 0.0_wp
           IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0_wp
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
-               precipitation )  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              sums_wsqrs_ws_l = 0.0_wp
              sums_wsnrs_ws_l = 0.0_wp

palm/trunk/SOURCE/boundary_conds.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-!
+! icloud_scheme removed. microphyisics_seifert added.
 !
 ! Former revisions:
@@ -139 +139 @@
                cloud_physics, dt_3d, humidity,                                 &
                ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_sa_t, ibc_uv_b,       &
-               ibc_uv_t, icloud_scheme, inflow_l, inflow_n, inflow_r, inflow_s,&
-               intermediate_timestep_count, large_scale_forcing, nest_domain,  &
-               nest_bound_l, nest_bound_s, nudging, ocean,                     &
-               outflow_l, outflow_n, outflow_r, outflow_s, passive_scalar,     &
-               precipitation, tsc, use_cmax
+               ibc_uv_t, inflow_l, inflow_n, inflow_r, inflow_s,               &
+               intermediate_timestep_count, large_scale_forcing,               &
+               microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s,  &
+               nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s,     &
+               passive_scalar, tsc, use_cmax
 
     USE grid_variables,                                                        &
@@ -319 +319 @@
        ENDIF
 
-       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 !             
 !--       Surface conditions rain water (Dirichlet)
@@ -373 +373 @@
        IF ( humidity  .OR.  passive_scalar )  THEN
           q_p(:,nys-1,:) = q_p(:,nys,:)
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-               precipitation)  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p(:,nys-1,:) = qr_p(:,nys,:)
              nr_p(:,nys-1,:) = nr_p(:,nys,:)
@@ -384 +383 @@
        IF ( humidity  .OR.  passive_scalar )  THEN
           q_p(:,nyn+1,:) = q_p(:,nyn,:)
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-               precipitation )  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
              nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
@@ -395 +393 @@
        IF ( humidity  .OR.  passive_scalar )  THEN
           q_p(:,:,nxl-1) = q_p(:,:,nxl)
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-               precipitation )  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
              nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
@@ -406 +403 @@
        IF ( humidity .OR. passive_scalar )  THEN
           q_p(:,:,nxr+1) = q_p(:,:,nxr)
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
              nr_p(:,:,nxr+1) = nr_p(:,:,nxr)

palm/trunk/SOURCE/calc_liquid_water_content.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! icloud_scheme removed. microphysics_seifert added.
 ! 
 ! Former revisions:
@@ -73 +73 @@
 
     USE control_parameters,                                                    &
-        ONLY:  icloud_scheme, precipitation
+        ONLY:  microphysics_seifert
 
     USE indices,                                                               &
@@ -122 +122 @@
 !
 !--          Compute the liquid water content
-             IF ( icloud_scheme == 0  .AND.  precipitation)  THEN
+             IF ( microphysics_seifert )  THEN
                 IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp ) THEN
                    qc(k,j,i) = q(k,j,i) - q_s - qr(k,j,i)
@@ -131 +131 @@
                    ql(k,j,i) = qr(k,j,i)
                 ENDIF
-             ELSEIF ( icloud_scheme == 0  .AND.  .NOT. precipitation )  THEN
+             ELSE
                 IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN
                    qc(k,j,i) = q(k,j,i) - q_s
                    ql(k,j,i) = qc(k,j,i)
                 ELSE
-                   qc(k,j,i) = 0.0_wp 
+                   qc(k,j,i) = 0.0_wp
                    ql(k,j,i) = 0.0_wp
-                ENDIF
-             ELSE
-                IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN
-                   ql(k,j,i) = q(k,j,i) - q_s
-                ELSE
-                   ql(k,j,i) = 0.0_wp 
                 ENDIF
              ENDIF

palm/trunk/SOURCE/check_parameters.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! PALM collision kernel deleted. Collision algorithms are checked for correct
+! spelling.
+!
+! Tails removed.
+!
+! Checks for use_sgs_for_particles adopted for the use of droplets with
+! use_sgs_for_particles adopted.
+!
+! Unused variables removed.
 ! 
 ! Former revisions:
@@ -358 +366 @@
     CHARACTER (LEN=1)   ::  sq                       !< 
     CHARACTER (LEN=15)  ::  var                      !< 
-    CHARACTER (LEN=7)   ::  unit                     !< 
+    CHARACTER (LEN=7)   ::  unit                     !<
     CHARACTER (LEN=8)   ::  date                     !< 
     CHARACTER (LEN=10)  ::  time                     !< 
@@ -367 +375 @@
 
     INTEGER(iwp) ::  i                               !< 
-    INTEGER(iwp) ::  ilen                            !< 
-    INTEGER(iwp) ::  iremote = 0                     !< 
+    INTEGER(iwp) ::  ilen                            !<
     INTEGER(iwp) ::  j                               !< 
     INTEGER(iwp) ::  k                               !< 
@@ -374 +381 @@
     INTEGER(iwp) ::  netcdf_data_format_save         !< 
     INTEGER(iwp) ::  position                        !<
-    INTEGER(iwp) ::  prec                            !< 
     
     LOGICAL     ::  found                            !< 
@@ -661 +667 @@
                         'with particle advection.'
        CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
-    ENDIF
-
-!
-!--
-    IF ( use_particle_tails )  THEN
-       message_string = 'Particle tails are currently not available due ' //   &
-                        'to the new particle structure.'
-       CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
     ENDIF
 
@@ -767 +765 @@
 !
 !-- Check cloud scheme
-    IF ( cloud_scheme == 'seifert_beheng' )  THEN
-       icloud_scheme = 0
+    IF ( cloud_scheme == 'saturation_adjust' )  THEN
+       microphysics_sat_adjust = .TRUE.
+       microphysics_seifert    = .FALSE.
+       microphysics_kessler    = .FALSE.
+       precipitation           = .FALSE.
+    ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
+       microphysics_sat_adjust = .FALSE.
+       microphysics_seifert    = .TRUE.
+       microphysics_kessler    = .FALSE.
+       precipitation           = .TRUE.
     ELSEIF ( cloud_scheme == 'kessler' )  THEN
-       icloud_scheme = 1
+       microphysics_sat_adjust = .FALSE.
+       microphysics_seifert    = .FALSE.
+       microphysics_kessler    = .TRUE.
+       precipitation           = .TRUE.
     ELSE
        message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
@@ -846 +855 @@
     ENDIF
 
-    IF ( use_sgs_for_particles  .AND.  .NOT. use_upstream_for_tke .AND.        &
-         ( scalar_advec /= 'ws-scheme' .OR.                                    &
-           scalar_advec /= 'ws-scheme-mono' )                                  &
+    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.               &
+         .NOT. use_upstream_for_tke  .AND.                                       &
+         ( scalar_advec /= 'ws-scheme'  .OR.  scalar_advec /= 'ws-scheme-mono' ) &
        )  THEN
        use_upstream_for_tke = .TRUE.
@@ -906 +915 @@
           hall_kernel = .TRUE.
 
-       CASE ( 'palm' )
-          palm_kernel = .TRUE.
-
        CASE ( 'wang', 'wang_fast' )
           wang_kernel = .TRUE.
@@ -922 +928 @@
     END SELECT
     IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
+
+!
+!-- Collision algorithms:
+    SELECT CASE ( TRIM( collision_algorithm ) )
+
+       CASE ( 'all_or_nothing' )
+          all_or_nothing = .TRUE.
+
+       CASE ( 'average_impact' )
+          average_impact = .TRUE.
+
+       CASE DEFAULT
+          message_string = 'unknown collision algorithm: collision_algorithm = "' // &
+                           TRIM( collision_algorithm ) // '"'
+          CALL message( 'check_parameters', 'PA0420', 1, 2, 0, 6, 0 )
+
+       END SELECT
 
     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
@@ -979 +1002 @@
               'not allowed with humidity = ', humidity
        CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
-    ENDIF
-
-    IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
-       WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ',   &
-              'not allowed with cloud_physics = ', cloud_physics
-       CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
     ENDIF
 
@@ -1032 +1049 @@
        ENDIF 
 
-       IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
           message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
                            ' seifert_beheng'
@@ -1095 +1112 @@
     IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
                  loop_optimization == 'vector' )                               &
-         .AND.  cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+         .AND.  cloud_physics  .AND.  microphysics_seifert )  THEN
        message_string = 'cloud_scheme = seifert_beheng requires ' //           &
                         'loop_optimization = "cache" or "vector"'
@@ -2152 +2169 @@
 !
 !-- Set the default value for the integration interval of precipitation amount
-    IF ( precipitation )  THEN
+    IF ( microphysics_seifert  .OR.  microphysics_kessler )  THEN
        IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
           precipitation_amount_interval = dt_do2d_xy
@@ -2782 +2799 @@
                                  'lemented for cloud_physics = .FALSE.'
                 CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
-             ELSEIF ( icloud_scheme /= 0 )  THEN
+             ELSEIF ( .NOT.  microphysics_seifert )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for cloud_scheme /= seifert_beheng'
                 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
-             ELSEIF (  .NOT.  precipitation )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for precipitation = .FALSE.'
-                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
              ELSE
                 dopr_index(i) = 73
@@ -2804 +2816 @@
                                  'lemented for cloud_physics = .FALSE.'
                 CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
-             ELSEIF ( icloud_scheme /= 0 )  THEN
+             ELSEIF ( .NOT.  microphysics_seifert )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for cloud_scheme /= seifert_beheng'
                 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
-             ELSEIF (  .NOT.  precipitation )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for precipitation = .FALSE.'
-                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
              ELSE
                 dopr_index(i) = 74
@@ -2826 +2833 @@
                                  'lemented for cloud_physics = .FALSE.'
                 CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
-             ELSEIF ( icloud_scheme /= 0 )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for cloud_scheme /= seifert_beheng'
-                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
              ELSE
                 dopr_index(i) = 75
@@ -2843 +2845 @@
                                  'lemented for cloud_physics = .FALSE.'
                 CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
-             ELSEIF ( icloud_scheme /= 0 )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for cloud_scheme /= seifert_beheng'
-                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
-             ELSEIF (  .NOT.  precipitation )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for precipitation = .FALSE.'
-                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
-
+             ELSEIF ( microphysics_sat_adjust )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'ilable for cloud_scheme = saturation_adjust'
+                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
              ELSE
                 dopr_index(i) = 76
@@ -3203 +3199 @@
                          'res cloud_physics = .TRUE.'
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
-             ELSEIF ( icloud_scheme /= 0 )  THEN
+             ELSEIF ( .NOT.  microphysics_seifert )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_scheme = seifert_beheng'
@@ -3224 +3220 @@
                          'res cloud_physics = .TRUE.'
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
-             ELSEIF ( icloud_scheme /= 0 )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                         'res cloud_scheme = seifert_beheng'
-                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
-             ELSEIF (  .NOT.  precipitation )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res precipitation = .TRUE.'
-                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
+             ELSEIF ( microphysics_sat_adjust )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" is ' //  &
+                         'not available for cloud_scheme = saturation_adjust'
+                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
              ENDIF
              unit = 'kg/kg m/s'
@@ -3249 +3241 @@
                          'res cloud_physics = .TRUE.'
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
-             ELSEIF ( icloud_scheme /= 0 ) THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                         'res cloud_scheme = seifert_beheng'
-                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
              ENDIF
              unit = 'kg/kg'
@@ -3279 +3267 @@
                          'res cloud_physics = .TRUE.'
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
-             ELSEIF ( icloud_scheme /= 0 ) THEN
+             ELSEIF ( .NOT.  microphysics_seifert ) THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_scheme = seifert_beheng'
                 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
-             ELSEIF (  .NOT.  precipitation )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res precipitation = .TRUE.'
-                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
              ENDIF
              unit = 'kg/kg'
@@ -3371 +3355 @@
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
+             IF ( TRIM( var ) == 'pra*'  .AND.                                 &
+                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res precipitation = .TRUE.'
+                                 'res cloud_scheme = kessler or seifert_beheng'
                 CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
              ENDIF
@@ -3381 +3366 @@
                 CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT.  precipitation )  THEN
+             IF ( TRIM( var ) == 'prr*'  .AND.                                 &
+                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res precipitation = .TRUE.'
+                                 'res cloud_scheme = kessler or seifert_beheng'
                 CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
              ENDIF

palm/trunk/SOURCE/data_output_2d.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of bulk cloud physics simplified.
 ! 
 ! Former revisions:
@@ -146 +146 @@
                do2d_xz_last_time, do2d_xz_n, do2d_xz_time_count,               &
                do2d_yz_last_time, do2d_yz_n, do2d_yz_time_count,               &
-               ibc_uv_b, icloud_scheme, io_blocks, io_group,                   &
-               message_string,                             &
-               ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz, psolver, section,        &
-               simulated_time,  simulated_time_chr, time_since_reference_point
+               ibc_uv_b, io_blocks, io_group, message_string,                  &
+               ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz,                          &
+               psolver, section, simulated_time, simulated_time_chr,           &
+               time_since_reference_point
         
     USE cpulog,                                                                &
@@ -683 +683 @@
 
              CASE ( 'prr*_xy' )        ! 2d-array
-                IF ( icloud_scheme == 1 )  THEN
-                   IF ( av == 0 )  THEN
-                      CALL exchange_horiz_2d( precipitation_rate )
-                      DO  i = nxlg, nxrg
-                         DO  j = nysg, nyng
-                            local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
+                IF ( av == 0 )  THEN
+                   CALL exchange_horiz_2d( prr(nzb+1,:,:) )
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = prr(nzb+1,j,i) * hyrho(nzb+1)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   CALL exchange_horiz_2d( prr_av(nzb+1,:,:) )
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = prr_av(nzb+1,j,i) * hyrho(nzb+1)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
+                IF ( av == 0 )  THEN
+                   CALL exchange_horiz( prr, nbgp )
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            local_pf(i,j,k) = prr(k,j,i) * hyrho(nzb+1)
                          ENDDO
                       ENDDO
-                   ELSE
-                      CALL exchange_horiz_2d( precipitation_rate_av )
-                      DO  i = nxlg, nxrg
-                         DO  j = nysg, nyng
-                            local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
-                         ENDDO
-                      ENDDO
-                   ENDIF
-                ELSE
-                   IF ( av == 0 )  THEN
-                      CALL exchange_horiz_2d( prr(nzb+1,:,:) )
-                      DO  i = nxlg, nxrg
-                         DO  j = nysg, nyng
-                            local_pf(i,j,nzb+1) = prr(nzb+1,j,i) * hyrho(nzb+1)
-                         ENDDO
-                      ENDDO
-                   ELSE
-                      CALL exchange_horiz_2d( prr_av(nzb+1,:,:) )
-                      DO  i = nxlg, nxrg
-                         DO  j = nysg, nyng
-                            local_pf(i,j,nzb+1) = prr_av(nzb+1,j,i) *          &
-                                                  hyrho(nzb+1)
-                         ENDDO
-                      ENDDO
-                   ENDIF
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
-                IF ( av == 0 )  THEN
-                   CALL exchange_horiz( prr, nbgp )
+                   ENDDO
+                ELSE
+                   CALL exchange_horiz( prr_av, nbgp )
                    DO  i = nxlg, nxrg
                       DO  j = nysg, nyng
                          DO  k = nzb, nzt+1
-                            local_pf(i,j,k) = prr(k,j,i)
-                         ENDDO
-                      ENDDO
-                   ENDDO
-                ELSE
-                   CALL exchange_horiz( prr_av, nbgp )
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         DO  k = nzb, nzt+1
-                            local_pf(i,j,k) = prr_av(k,j,i)
+                            local_pf(i,j,k) = prr_av(k,j,i) * hyrho(nzb+1)
                          ENDDO
                       ENDDO

palm/trunk/SOURCE/data_output_3d.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-!
+! prr vertical dimensions set to nzb_do to nzt_do. Unused variables deleted.
 !
 ! Former revisions:
@@ -122 +122 @@
         
     USE control_parameters,                                                    &
-        ONLY:  avs_data_file, cloud_physics, do3d, do3d_avs_n,                 &
-               do3d_no, do3d_time_count, io_blocks, io_group,                  &
-               message_string, ntdim_3d,                                       &
-               nz_do3d, plot_3d_precision, psolver, simulated_time,            &
-               simulated_time_chr, skip_do_avs, time_since_reference_point
+        ONLY:  cloud_physics, do3d, do3d_no, do3d_time_count, io_blocks,       &
+               io_group, message_string, ntdim_3d, nz_do3d, psolver,           &
+               simulated_time, time_since_reference_point
         
     USE cpulog,                                                                &
@@ -132 +130 @@
         
     USE indices,                                                               &
-        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzt,  &
-               nzb
+        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb
         
     USE kinds
@@ -163 +160 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=9) ::  simulated_time_mod  !< 
-
     INTEGER(iwp) ::  av        !< 
     INTEGER(iwp) ::  i         !< 
@@ -173 +168 @@
     INTEGER(iwp) ::  nzb_do    !< vertical lower limit for data output
     INTEGER(iwp) ::  nzt_do    !< vertical upper limit for data output
-    INTEGER(iwp) ::  pos       !< 
-    INTEGER(iwp) ::  prec      !< 
-    INTEGER(iwp) ::  psi       !< 
 
     LOGICAL      ::  found     !< 
@@ -374 +366 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          local_pf(i,j,k) = prr(k,j,i)
                       ENDDO
@@ -383 +375 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          local_pf(i,j,k) = prr_av(k,j,i)
                       ENDDO

palm/trunk/SOURCE/data_output_dvrp.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Particles and tails removed.
 ! 
 ! Former revisions:
@@ -58 +58 @@
 ! Description:
 ! ------------
-!> Plot of isosurface, particles and slicers with dvrp-software
+!> Plot of isosurface and slicers with dvrp-software
 !------------------------------------------------------------------------------!
  MODULE dvrp_color
@@ -130 +130 @@
         
     USE kinds
-    
-    USE particle_attributes,                                                   &
-        ONLY:  maximum_number_of_tailpoints, number_of_particles,              &
-               number_of_tails, particle_advection, particle_advection_start,  &
-               particle_tail_coordinates, particles, uniform_particles,        &
-               use_particle_tails
-        
+
     USE pegrid
 
@@ -234 +228 @@
 !--    Select the plot mode (in case of isosurface or slicer only if user has
 !--    defined a variable which shall be plotted; otherwise do nothing)
-       IF ( mode_dvrp(m)(1:9) == 'particles'  .AND.  particle_advection  .AND. &
-            simulated_time >= particle_advection_start )  THEN
-!
-!--       DVRP-Calls for plotting particles:
-          CALL cpu_log( log_point_s(28), 'dvrp_particles', 'start' )
-
-!
-!--       Definition of characteristics of particle material
-!          tmp_r = 0.1;  tmp_g = 0.7;  tmp_b = 0.1;  tmp_t = 0.0
-          tmp_r = 0.0_wp;  tmp_g = 0.0_wp;  tmp_b = 0.0_wp;  tmp_t = 0.0_wp
-          CALL DVRP_MATERIAL_RGB( m-1, 1, tmp_r, tmp_g, tmp_b, tmp_t )
-
-!
-!--       If clipping is active and if this subdomain is clipped, find out the
-!--       number of particles and tails to be plotted; otherwise, all
-!--       particles/tails are plotted. dvrp_mask is used to mark the partikles.
-          ALLOCATE( dvrp_mask(number_of_particles) )
-
-          IF ( dvrp_total_overlap )  THEN
-             dvrp_mask = .TRUE.
-             dvrp_nop  = number_of_particles
-             dvrp_not  = number_of_tails
-          ELSE
-             dvrp_mask = .FALSE.
-             dvrp_nop  = 0
-             dvrp_not  = 0
-             IF ( dvrp_overlap )  THEN
-                IF ( .NOT. use_particle_tails )  THEN
-                   DO  n = 1, number_of_particles
-                      ip = particles(n)%x / dx
-                      jp = particles(n)%y / dy
-                      IF ( ip >= nxl_dvrp  .AND.  ip <= nxr_dvrp  .AND.  &
-                           jp >= nys_dvrp  .AND.  jp <= nyn_dvrp )  THEN
-                         dvrp_nop = dvrp_nop + 1
-                         dvrp_mask(n) = .TRUE.
-                      ENDIF
-                   ENDDO
-                ELSE
-                   k = 0
-                   DO  n = 1, number_of_particles
-                      IF ( particles(n)%tail_id /= 0 )  THEN
-                         k = k + 1
-                         ip = particles(n)%x / dx
-                         jp = particles(n)%y / dy
-                         IF ( ip >= nxl_dvrp  .AND.  ip <= nxr_dvrp  .AND.  &
-                              jp >= nys_dvrp  .AND.  jp <= nyn_dvrp )  THEN
-                            dvrp_not = dvrp_not + 1
-                            dvrp_mask(n) = .TRUE.
-                         ENDIF
-                      ENDIF
-                   ENDDO
-                ENDIF
-             ENDIF
-          ENDIF
-
-!
-!--       Move particle coordinates to one-dimensional arrays
-          IF ( .NOT. use_particle_tails )  THEN
-!
-!--          All particles are output
-             ALLOCATE( psize(dvrp_nop), p_t(dvrp_nop), p_c(dvrp_nop), &
-                       p_x(dvrp_nop), p_y(dvrp_nop), p_z(dvrp_nop) )
-             psize = 0.0_wp;  p_t = 0_wp;  p_c = 0.0_wp
-             p_x = 0.0_wp; p_y = 0.0_wp
-             p_z   = 0.0_wp
-             k = 0
-             DO  n = 1, number_of_particles
-                IF ( dvrp_mask(n) )  THEN
-                   k = k + 1
-                   psize(k) = particles(n)%dvrp_psize
-                   p_x(k)   = particles(n)%x * superelevation_x
-                   p_y(k)   = particles(n)%y * superelevation_y
-                   p_z(k)   = particles(n)%z * superelevation
-                   p_c(k)   = particles(n)%class
-                ENDIF
-             ENDDO
-          ELSE
-!
-!--          Particles have a tail
-             ALLOCATE( psize(dvrp_not), p_t(dvrp_not),             &
-                       p_c(dvrp_not*maximum_number_of_tailpoints), &
-                       p_x(dvrp_not*maximum_number_of_tailpoints), &
-                       p_y(dvrp_not*maximum_number_of_tailpoints), &
-                       p_z(dvrp_not*maximum_number_of_tailpoints) )
-             psize = 0.0_wp;  p_t = 0_wp;  p_c = 0.0_wp
-             p_x = 0.0_wp;  p_y = 0.0_wp
-             p_z   = 0.0_wp
-             i = 0
-             k = 0
-
-             DO  n = 1, number_of_particles
-                nn = particles(n)%tail_id
-                IF ( nn /= 0  .AND.  dvrp_mask(n) )  THEN
-                   k = k + 1
-                   DO  j = 1, particles(n)%tailpoints
-                      i = i + 1
-                      p_x(i) = particle_tail_coordinates(j,1,nn) * &
-                                                                superelevation_x
-                      p_y(i) = particle_tail_coordinates(j,2,nn) * &
-                                                                superelevation_y
-                      p_z(i) = particle_tail_coordinates(j,3,nn) * &
-                                                                superelevation
-                      p_c(i) = particle_tail_coordinates(j,4,nn)
-                   ENDDO
-                   psize(k) = particles(n)%dvrp_psize
-                   p_t(k)   = particles(n)%tailpoints - 1
-                ENDIF                
-             ENDDO
-          ENDIF
-
-!
-!--       Compute and plot particles in dvr-format
-          IF ( uniform_particles  .AND.  .NOT. use_particle_tails )  THEN
-!
-!--          All particles have the same color. Use simple routine to set
-!--          the particle attributes (produces less output data)
-             CALL DVRP_PARTICLES( m-1, p_x, p_y, p_z, psize )
-          ELSE
-!
-!--          Set color definitions
-             CALL user_dvrp_coltab( 'particles', 'none' )
-
-             CALL DVRP_COLORTABLE_HLS( m-1, 0, interval_values_dvrp_prt, &
-                                       interval_h_dvrp_prt,              &
-                                       interval_l_dvrp_prt,              &
-                                       interval_s_dvrp_prt,              &
-                                       interval_a_dvrp_prt )
-
-             IF ( .NOT. use_particle_tails )  THEN
-                CALL DVRP_PARTICLES( m-1, dvrp_nop, p_x, p_y, p_z, 3, psize, &
-                                     p_c, p_t )
-             ELSE
-                CALL DVRP_PARTICLES( m-1, dvrp_not, p_x, p_y, p_z, 15, psize, &
-                                     p_c, p_t )
-             ENDIF
-          ENDIF
-
-          CALL DVRP_VISUALIZE( m-1, 3, dvrp_filecount )
-
-          DEALLOCATE( dvrp_mask, psize, p_c, p_t, p_x, p_y, p_z )
-
-          CALL cpu_log( log_point_s(28), 'dvrp_particles', 'stop' )
-
-
-       ELSEIF ( ( mode_dvrp(m)(1:10) == 'isosurface'  .OR.                     &
+       IF ( ( mode_dvrp(m)(1:10) == 'isosurface'  .OR.                     &
                   mode_dvrp(m)(1:6)  == 'slicer'           )                   &
                   .AND.  output_variable /= ' ' )  THEN

palm/trunk/SOURCE/flow_statistics.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-!
+! Output of bulk microphysics simplified.
 !
 ! Former revisions:
@@ -193 +193 @@
     USE control_parameters,                                                    &
         ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
-                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
+                dt_3d, g, humidity, kappa, large_scale_forcing,                &
                 large_scale_subsidence, max_pr_user, message_string, neutral,  &
-                ocean, passive_scalar, precipitation, simulated_time,          &
+                microphysics_seifert, ocean, passive_scalar, simulated_time,   &
                 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
                 ws_scheme_mom, ws_scheme_sca
@@ -898 +898 @@
                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
                                                           rmask(j,i,sr)
+                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
+
                       IF ( .NOT. cloud_droplets )  THEN
                          pts = 0.5_wp *                                        &
@@ -906 +908 @@
                          sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
                                                              rmask(j,i,sr)
-                         IF ( icloud_scheme == 0  )  THEN
-                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
+                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
+                                                             rmask(j,i,sr)
+                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
+                                                             rmask(j,i,sr)
+                         IF ( microphysics_seifert )  THEN
+                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
                                                                 rmask(j,i,sr)
-                            sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *    &
-                                                                rmask(j,i,sr)
-                            IF ( precipitation )  THEN
-                               sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * &
-                                                                   rmask(j,i,sr)
-                               sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * &
-                                                                   rmask(j,i,sr)
-                               sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *&
-                                                                   rmask(j,i,sr)
-                            ENDIF
-                         ELSE
-                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
+                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
                                                                 rmask(j,i,sr)
                          ENDIF
-                      ELSE
-                         sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *       &
-                                                             rmask(j,i,sr)
                       ENDIF
+
                    ELSE
                       IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
@@ -1662 +1655 @@
     USE control_parameters,                                                    &
         ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
-                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
-                large_scale_subsidence, max_pr_user, message_string, neutral,  &
-                ocean, passive_scalar, precipitation, simulated_time,          &
-                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
-                ws_scheme_mom, ws_scheme_sca
+                dt_3d, g, humidity, kappa, large_scale_forcing,                &
+                large_scale_subsidence, max_pr_user, message_string,           &
+                microphysics_seifert, neutral, ocean, passive_scalar,          &
+                simulated_time, use_subsidence_tendencies, use_surface_fluxes, &
+                use_top_fluxes, ws_scheme_mom, ws_scheme_sca
         
     USE cpulog,                                                                &
@@ -2774 +2767 @@
                 !$acc end parallel loop
 
-                IF ( icloud_scheme == 0  )  THEN
+                IF ( microphysics_seifert )  THEN
 
                    !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
@@ -2791 +2784 @@
                    !$acc end parallel loop
 
-                   IF ( precipitation )  THEN
-
-                      !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
-                      DO  k = nzb, nzt_diff
-                         s1 = 0
-                         s2 = 0
-                         s3 = 0
-                         !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
-                         DO  i = nxl, nxr
-                            DO  j = nys, nyn
-                               s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
-                               s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
-                               s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
-                            ENDDO
+                   !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
+                   DO  k = nzb, nzt_diff
+                      s1 = 0
+                      s2 = 0
+                      s3 = 0
+                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
+                      DO  i = nxl, nxr
+                         DO  j = nys, nyn
+                            s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                            s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                            s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                          ENDDO
-                         sums_l(k,73,tn) = s1
-                         sums_l(k,74,tn) = s2
-                         sums_l(k,76,tn) = s3
                       ENDDO
-                      !$acc end parallel loop
-
-                   ENDIF
+                      sums_l(k,73,tn) = s1
+                      sums_l(k,74,tn) = s2
+                      sums_l(k,76,tn) = s3
+                   ENDDO
+                   !$acc end parallel loop
 
                 ELSE

palm/trunk/SOURCE/header.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+! Tails removed. icloud_scheme replaced by microphysics_*
+!
 ! Former revisions:
 ! -----------------
@@ -275 +275 @@
                density_ratio, dissipation_classes, dt_min_part, dt_prel,       &
                dt_write_particle_data, end_time_prel,                          &
-               maximum_number_of_tailpoints, maximum_tailpoint_age,            &
-               minimum_tailpoint_distance, number_of_particle_groups,          &
-               particle_advection, particle_advection_start,                   &
+               number_of_particle_groups, particle_advection,                  &
+               particle_advection_start,                                       &
                particles_per_point, pdx, pdy, pdz,  psb, psl, psn, psr, pss,   &
                pst, radius, radius_classes, random_start_position,             &
                seed_follows_topography,                                        &
-               total_number_of_particles, use_particle_tails,                  &
-               use_sgs_for_particles, total_number_of_tails,                   &
+               total_number_of_particles, use_sgs_for_particles,               &
                vertical_particle_advection, write_particle_statistics
         
@@ -1700 +1698 @@
                                    slicer_range_limits_dvrp(:,m)
              ENDIF
-          ELSEIF ( mode_dvrp(i)(1:9) == 'particles' )  THEN
-             WRITE ( io, 363 )  dvrp_psize
-             IF ( particle_dvrpsize /= 'none' )  THEN
-                WRITE ( io, 364 )  'size', TRIM( particle_dvrpsize ), &
-                                   dvrpsize_interval
-             ENDIF
-             IF ( particle_color /= 'none' )  THEN
-                WRITE ( io, 364 )  'color', TRIM( particle_color ), &
-                                   color_interval
-             ENDIF
           ENDIF
           i = i + 1
@@ -1782 +1770 @@
        WRITE ( io, 415 )
        WRITE ( io, 416 ) surface_pressure, r_d, rho_surface, cp, l_v
-       IF ( icloud_scheme == 0 )  THEN
+       IF ( microphysics_seifert )  THEN
           WRITE ( io, 510 ) 1.0E-6_wp * nc_const
-          IF ( precipitation )  WRITE ( io, 511 ) c_sedimentation
+          WRITE ( io, 511 ) c_sedimentation
        ENDIF
     ENDIF
@@ -1796 +1784 @@
        WRITE ( io, 432 )
        IF ( cloud_top_radiation )  WRITE ( io, 132 )
-       IF ( icloud_scheme == 1 )  THEN
-          IF ( precipitation )  WRITE ( io, 133 )
-       ELSEIF ( icloud_scheme == 0 )  THEN
+       IF ( microphysics_kessler )  THEN
+          WRITE ( io, 133 )
+       ELSEIF ( microphysics_seifert )  THEN
           IF ( drizzle )  WRITE ( io, 506 )
-          IF ( precipitation )  THEN
-             WRITE ( io, 505 )
-             IF ( turbulence )  WRITE ( io, 507 )
-             IF ( ventilation_effect )  WRITE ( io, 508 )
-             IF ( limiter_sedimentation )  WRITE ( io, 509 )
-          ENDIF
+          WRITE ( io, 505 )
+          IF ( turbulence )  WRITE ( io, 507 )
+          IF ( ventilation_effect )  WRITE ( io, 508 )
+          IF ( limiter_sedimentation )  WRITE ( io, 509 )
        ENDIF
     ELSEIF ( humidity  .AND.  cloud_droplets )  THEN
@@ -1872 +1858 @@
        IF ( particles_per_point > 1 )  WRITE ( io, 489 )  particles_per_point
        WRITE ( io, 495 )  total_number_of_particles
-       IF ( use_particle_tails  .AND.  maximum_number_of_tailpoints /= 0 )  THEN
-          WRITE ( io, 483 )  maximum_number_of_tailpoints
-          IF ( minimum_tailpoint_distance /= 0 )  THEN
-             WRITE ( io, 484 )  total_number_of_tails,      &
-                                minimum_tailpoint_distance, &
-                                maximum_tailpoint_age
-          ENDIF
-       ENDIF
        IF ( dt_write_particle_data /= 9999999.9_wp )  THEN
           WRITE ( io, 485 )  dt_write_particle_data
@@ -2215 +2193 @@
 362 FORMAT (/'       Slicer plane ',A/ &
             '       Slicer limits: [',F6.2,',',F6.2,']')
-363 FORMAT (/'       Particles'/ &
-            '          particle size:  ',F7.2,' m')
-364 FORMAT ('          particle ',A,' controlled by "',A,'" with interval [', &
-                       F6.2,',',F6.2,']')
 365 FORMAT (/'       Groundplate color: (',F4.2,',',F4.2,',',F4.2,') (R,G,B)'/ &
             '       Superelevation along (x,y,z): (',F4.1,',',F4.1,',',F4.1, &
@@ -2377 +2351 @@
 481 FORMAT ('       Particles have random start positions'/)
 482 FORMAT ('          Particles are advected only horizontally'/)
-483 FORMAT ('       Particles have tails with a maximum of ',I3,' points')
-484 FORMAT ('            Number of tails of the total domain: ',I10/ &
-            '            Minimum distance between tailpoints: ',F8.2,' m'/ &
-            '            Maximum age of the end of the tail:  ',F8.2,' s')
 485 FORMAT ('       Particle data are written on file every ', F9.1, ' s')
 486 FORMAT ('       Particle statistics are written on file'/)

palm/trunk/SOURCE/init_3d_model.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! icloud_scheme replaced by microphysics_*
 ! 
 ! Former revisions:
@@ -465 +465 @@
                        precipitation_rate(nysg:nyng,nxlg:nxrg) )
 
-             IF ( icloud_scheme == 0 )  THEN
-!
-!--             1D-arrays 
-                ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1),                 &
-                           q_1d(nzb:nzt+1), qc_1d(nzb:nzt+1) ) 
-!
-!--             3D-cloud water content
+!
+!--          1D-arrays
+             ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1),                 &
+                        q_1d(nzb:nzt+1), qc_1d(nzb:nzt+1) )
+!
+!--          3D-cloud water content
 #if defined( __nopointer )
-                ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #else
-                ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #endif
-
-                IF ( precipitation )  THEN
-!
-!--                1D-arrays 
-                   ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) ) 
-
-!
-!--                2D-rain water content and rain drop concentration arrays
-                   ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                        &
-                              qrsws(nysg:nyng,nxlg:nxrg),                      &
-                              qrswst(nysg:nyng,nxlg:nxrg),                     &
-                              nrs(nysg:nyng,nxlg:nxrg),                        &
-                              nrsws(nysg:nyng,nxlg:nxrg),                      &
-                              nrswst(nysg:nyng,nxlg:nxrg) )
-!
-!--                3D-rain water content, rain drop concentration arrays
+!
+!--          3d-precipitation rate
+             ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+             IF ( microphysics_seifert )  THEN
+!
+!--             1D-arrays
+                ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
+
+!
+!--             2D-rain water content and rain drop concentration arrays
+                ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                        &
+                           qrsws(nysg:nyng,nxlg:nxrg),                      &
+                           qrswst(nysg:nyng,nxlg:nxrg),                     &
+                           nrs(nysg:nyng,nxlg:nxrg),                        &
+                           nrsws(nysg:nyng,nxlg:nxrg),                      &
+                           nrswst(nysg:nyng,nxlg:nxrg) )
+!
+!--             3D-rain water content, rain drop concentration arrays
 #if defined( __nopointer )
-                   ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
-                             nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
-                             qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
-                             qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
-                             tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             &
-                             tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
+                          nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
+                          qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             &
+                          tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #else
-                   ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
-                             nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
-                             nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
-                             qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
-                             qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
-                             qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #endif
-!
-!--                3d-precipitation rate
-                   ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-
-             ENDIF
+             ENDIF
+
           ENDIF
 
@@ -658 +656 @@
           IF ( cloud_physics )  THEN
              ql => ql_1
-             IF ( icloud_scheme == 0 )  THEN
-                qc => qc_1
-                IF ( precipitation )  THEN
-                   qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
-                   nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
-                ENDIF
+             qc => qc_1
+             IF ( microphysics_seifert )  THEN
+                qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
+                nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
              ENDIF
           ENDIF
@@ -733 +729 @@
                 ENDDO
              ENDDO
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-                  precipitation )  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -790 +785 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              qs = 0.0_wp
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-                  precipitation )  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 qrs = 0.0_wp
                 nrs = 0.0_wp
@@ -880 +874 @@
                 ENDDO
              ENDDO
-             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 !
 !--             Initialze nc_1d with default value
                 nc_1d(:) = nc_const
 
-                IF ( precipitation )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         qr(:,j,i) = 0.0_wp
-                         nr(:,j,i) = 0.0_wp
-                      ENDDO
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr(:,j,i) = 0.0_wp
+                      nr(:,j,i) = 0.0_wp
                    ENDDO
-                ENDIF
+                ENDDO
 
              ENDIF
@@ -1078 +1070 @@
 !--       Determine the near-surface water flux
           IF ( humidity  .OR.  passive_scalar )  THEN
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-                  precipitation )  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 qrsws = 0.0_wp
                 nrsws = 0.0_wp
@@ -1116 +1107 @@
              IF ( humidity  .OR.  passive_scalar )  THEN
                 qswst = 0.0_wp
-                IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.           &
-                     precipitation ) THEN
+                IF ( cloud_physics  .AND.  microphysics_seifert ) THEN
                    nrswst = 0.0_wp
                    qrswst = 0.0_wp
@@ -1157 +1147 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              IF (  .NOT.  constant_waterflux )  qsws   = 0.0_wp
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-                  precipitation )  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 qrsws = 0.0_wp
                 nrsws = 0.0_wp
@@ -1200 +1189 @@
        IF ( cloud_physics )  THEN
           ql = 0.0_wp
-          IF ( precipitation )  precipitation_amount = 0.0_wp
-          IF ( icloud_scheme == 0 )  THEN
-             qc = 0.0_wp
+          qc = 0.0_wp
+
+          precipitation_amount = 0.0_wp
+
+          IF ( microphysics_seifert )  THEN
              nc_1d = nc_const
           ENDIF 
@@ -1240 +1231 @@
           tq_m = 0.0_wp
           q_p = q
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
-               precipitation )  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              tqr_m = 0.0_wp
-             qr_p = qr
+             qr_p  = qr
              tnr_m = 0.0_wp
-             nr_p = nr
+             nr_p  = nr
           ENDIF
        ENDIF
@@ -1426 +1416 @@
        IF ( humidity  .OR.  passive_scalar )  THEN
           q_p = q
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
-               precipitation )  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p = qr
              nr_p = nr
@@ -1441 +1430 @@
        IF ( humidity  .OR.  passive_scalar )  THEN
           tq_m = 0.0_wp
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
-               precipitation )  THEN
+          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              tqr_m = 0.0_wp
              tnr_m = 0.0_wp

palm/trunk/SOURCE/init_cloud_physics.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! icloud_scheme replaced by microphysics_*
 ! 
 ! Former revisions:
@@ -92 +92 @@
         
     USE control_parameters,                                                    &
-        ONLY:  g, icloud_scheme, message_string, precipitation, pt_surface,    &
+        ONLY:  g, message_string, microphysics_seifert, pt_surface,            &
                rho_surface, surface_pressure
     
@@ -126 +126 @@
 !
 !-- Calculate timestep according to precipitation
-    IF ( icloud_scheme == 0  .AND.  precipitation )  THEN
+    IF ( microphysics_seifert )  THEN
        dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /         &
                           w_precipitation

palm/trunk/SOURCE/init_masks.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! icloud_scheme replaced by microphysics_*
 ! 
 ! Former revisions:
@@ -97 +97 @@
                data_output_masks, data_output_masks_user,                      &
                doav, doav_n, domask, domask_no, dz, dz_stretch_level, humidity,&
-               icloud_scheme, mask, masks, mask_scale, mask_i,                 &
+               mask, masks, mask_scale, mask_i,                                &
                mask_i_global, mask_j, mask_j_global, mask_k, mask_k_global,    &
                mask_loop, mask_size, mask_size_l, mask_start_l, mask_x,        &
                mask_x_loop, mask_xyz_dimension, mask_y, mask_y_loop, mask_z,   &
                mask_z_loop, max_masks,  message_string, mid,                   &
-               passive_scalar, precipitation, ocean
+               microphysics_seifert, passive_scalar, ocean
 
     USE grid_variables,                                                        &
@@ -270 +270 @@
                         '" requires cloud_physics = .TRUE.'
                    CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
-                 ELSEIF ( icloud_scheme /= 0 ) THEN
+                 ELSEIF ( .NOT. microphysics_seifert ) THEN
                    message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_scheme = seifert_beheng'
                    CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
-                 ELSEIF ( .NOT. precipitation )  THEN
-                   message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res precipitation = .TRUE.'
-                   CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
                 ENDIF
                 unit = '1/m3'
@@ -305 +301 @@
                             'res cloud_physics = .TRUE.'
                    CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
-                ELSEIF ( icloud_scheme /= 0 ) THEN
-                   message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                            'res cloud_scheme = seifert_beheng'
-                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
                 ENDIF
                 unit = 'kg/kg'
@@ -344 +336 @@
                             'res cloud_physics = .TRUE.'
                    CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
-                ELSEIF ( icloud_scheme /= 0 ) THEN
+                ELSEIF ( .NOT. microphysics_seifert ) THEN
                    message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                             'res cloud_scheme = seifert_beheng'
                    CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
-                ELSEIF ( .NOT. precipitation )  THEN
-                   message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                    'res precipitation = .TRUE.'
-                   CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
                 ENDIF
                 unit = 'kg/kg'

palm/trunk/SOURCE/interaction_droplets_ptq.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Unused variables removed.
 !
 ! Former revisions:
@@ -127 +127 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nyn, nys, nzb_2d, nzt
+           ONLY:  nzb_2d, nzt
 
        USE kinds,                                                              &

palm/trunk/SOURCE/lpm.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Tails removed.
+!
+! Initialization of sgs model not necessary for the use of cloud_droplets and
+! use_sgs_for_particles.
+!
+! lpm_release_set integrated.
+!
+! Unused variabled removed.
 ! 
 ! Former revisions:
@@ -108 +115 @@
         ONLY:  lpm_exchange_horiz, lpm_move_particle
 
+    USE lpm_init_mod,                                                          &
+        ONLY: lpm_create_particle, PHASE_RELEASE
+
     USE lpm_pack_arrays_mod,                                                   &
         ONLY:  lpm_pack_all_arrays
 
     USE particle_attributes,                                                   &
-        ONLY:  collision_kernel, deleted_particles, deleted_tails,             &
+        ONLY:  collision_kernel, deleted_particles,                            &
                dt_write_particle_data, dt_prel, end_time_prel,                 &
                grid_particles, number_of_particles, number_of_particle_groups, &
                particles, particle_groups, prt_count, trlp_count_sum,          &
-               tail_mask, time_prel, time_sort_particles,                      &
+               time_prel,                                                      &
                time_write_particle_data, trlp_count_recv_sum, trnp_count_sum,  &
                trnp_count_recv_sum, trrp_count_sum, trrp_count_recv_sum,       &
-               trsp_count_sum, trsp_count_recv_sum, use_particle_tails,        &
+               trsp_count_sum, trsp_count_recv_sum,                            &
                use_sgs_for_particles, write_particle_statistics
 
@@ -159 +169 @@
 
 !
-!-- Initialize arrays for marking those particles/tails to be deleted after the
+!-- Initialize arrays for marking those particles to be deleted after the
 !-- (sub-) timestep
     deleted_particles = 0
-
-    IF ( use_particle_tails )  THEN
-       tail_mask = .TRUE.
-    ENDIF
-    deleted_tails = 0
-
 
 !
@@ -202 +206 @@
     IF ( time_prel >= dt_prel  .AND.  end_time_prel > simulated_time )  THEN
 
-       CALL lpm_release_set
+       CALL lpm_create_particle(PHASE_RELEASE)
 !
 !--    The MOD function allows for changes in the output interval with
@@ -240 +244 @@
 !--    velocity variances)
 
-       IF ( use_sgs_for_particles )  CALL lpm_init_sgs_tke
+       IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
+          CALL lpm_init_sgs_tke
+       ENDIF
 
        DO  i = nxl, nxr
@@ -361 +367 @@
        deleted_particles = 0
 
-       IF ( use_particle_tails )  THEN
-          tail_mask     = .TRUE.
-       ENDIF
-       deleted_tails = 0
-
-
        IF ( dt_3d_reached )  EXIT
 
@@ -387 +387 @@
     IF ( collision_kernel == 'none' )  CALL lpm_set_attributes
 
-
 !
 !-- Set particle attributes defined by the user
     CALL user_lpm_set_attributes
-
-
-!
-!-- particle tails currently not available
-!
-!-- If required, add the current particle positions to the particle tails
-!   IF ( use_particle_tails )  CALL lpm_extend_tails
-
 
 !

palm/trunk/SOURCE/lpm_advec.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Random velocity fluctuations for particles added. Terminal fall velocity
+! for droplets is calculated from a parameterization (which is better than
+! the previous, physically correct calculation, which demands a very short
+! time step that is not used in the model).
+!
+! Unused variables deleted.
 ! 
 ! Former revisions:
@@ -78 +83 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  de_dx, de_dy, de_dz, diss, e, u, us, usws, v, vsws, w, z0, zu,  &
-               zw
+        ONLY:  de_dx, de_dy, de_dz, diss, e, km, u, us, usws, v, vsws, w, zu, zw
 
     USE cpulog
@@ -87 +91 @@
     USE control_parameters,                                                    &
         ONLY:  atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d,   &
-               dt_3d_reached_l, dz, g, kappa, molecular_viscosity, topography, &
-               u_gtrans, v_gtrans, simulated_time
+               dt_3d_reached_l, dz, g, kappa, topography, u_gtrans, v_gtrans
 
     USE grid_variables,                                                        &
@@ -99 +102 @@
     
     USE particle_attributes,                                                   &
-        ONLY:  block_offset, c_0, density_ratio, dt_min_part, grid_particles,  &
+        ONLY:  block_offset, c_0, dt_min_part, grid_particles,                 &
                iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
                particles, particle_groups, offset_ocean_nzt, sgs_wfu_part,     &
@@ -140 +143 @@
     REAL(wp) ::  de_dz_int_l        !<
     REAL(wp) ::  de_dz_int_u        !<
+    REAL(wp) ::  diameter           !< diamter of droplet
     REAL(wp) ::  diss_int_l         !<
     REAL(wp) ::  diss_int_u         !<
@@ -150 +154 @@
     REAL(wp) ::  exp_term           !<
     REAL(wp) ::  gg                 !<
-    REAL(wp) ::  height_int         !<
     REAL(wp) ::  height_p           !<
-    REAL(wp) ::  lagr_timescale     !<
+    REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
     REAL(wp) ::  location(1:30,1:3) !<
     REAL(wp) ::  random_gauss       !<
+    REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
+    REAL(wp) ::  rg1                !< Gaussian distributed random number
+    REAL(wp) ::  rg2                !< Gaussian distributed random number
+    REAL(wp) ::  rg3                !< Gaussian distributed random number
+    REAL(wp) ::  sigma              !< velocity standard deviation
     REAL(wp) ::  u_int_l            !<
     REAL(wp) ::  u_int_u            !<
@@ -163 +171 @@
     REAL(wp) ::  w_int_l            !<
     REAL(wp) ::  w_int_u            !<
+    REAL(wp) ::  w_s                !< terminal velocity of droplets
     REAL(wp) ::  x                  !<
     REAL(wp) ::  y                  !<
-    REAL(wp) ::  z_p                !<    
+    REAL(wp) ::  z_p                !<
+
+    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
+    REAL(wp), PARAMETER ::  b_rog = 10.43_wp     !< parameter for fall velocity
+    REAL(wp), PARAMETER ::  c_rog = 0.6_wp       !< parameter for fall velocity
+    REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp   !< parameter for fall velocity
+    REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp  !< parameter for fall velocity
+    REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
 
     REAL(wp), DIMENSION(1:30) ::  d_gp_pl !<
@@ -388 +404 @@
 !-- Interpolate and calculate quantities needed for calculating the SGS
 !-- velocities
-    IF ( use_sgs_for_particles )  THEN
+    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
 
        IF ( topography == 'flat' )  THEN
@@ -1336 +1352 @@
 !--          Update of the particle velocity
              IF ( cloud_droplets )  THEN
-                exp_arg  = 4.5_wp * dens_ratio(n) * molecular_viscosity /      &
-                           ( particles(n)%radius )**2 *                        &
-                           ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius &
-                             * SQRT( ( u_int(n) - particles(n)%speed_x )**2 +  &
-                                     ( v_int(n) - particles(n)%speed_y )**2 +  &
-                                     ( w_int(n) - particles(n)%speed_z )**2 )  &
-                                          / molecular_viscosity )**0.687_wp    &
-                           )
-
-                exp_term = EXP( -exp_arg * dt_particle(n) )
-             ELSEIF ( use_sgs_for_particles )  THEN
+!
+!--             Terminal velocity is computed for vertical direction (Rogers et 
+!--             al., 1993, J. Appl. Meteorol.)
+                diameter = particles(n)%radius * 2000.0_wp !diameter in mm
+                IF ( diameter <= d0_rog )  THEN
+                   w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
+                ELSE
+                   w_s = a_rog - b_rog * EXP( -c_rog * diameter )
+                ENDIF
+
+!
+!--             If selected, add random velocities following Soelch and Kaercher
+!--             (2010, Q. J. R. Meteorol. Soc.)
+                IF ( use_sgs_for_particles )  THEN
+                   lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
+                   RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
+                   sigma          = SQRT( e(kp,jp,ip) )
+
+                   rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+                   rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+                   rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+
+                   particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
+                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg1
+                   particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
+                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg2
+                   particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
+                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg3
+
+                   particles(n)%speed_x = u_int(n) + particles(n)%rvar1
+                   particles(n)%speed_y = v_int(n) + particles(n)%rvar2
+                   particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
+                ELSE
+                   particles(n)%speed_x = u_int(n)
+                   particles(n)%speed_y = v_int(n)
+                   particles(n)%speed_z = w_int(n) - w_s
+                ENDIF
+
+             ELSE
+
+                IF ( use_sgs_for_particles )  THEN
+                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
+                   exp_term = EXP( -exp_arg * dt_particle(n) )
+                ELSE
+                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
+                   exp_term = particle_groups(particles(n)%group)%exp_term
+                ENDIF
+                particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
+                                       u_int(n) * ( 1.0_wp - exp_term )
+                particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
+                                       v_int(n) * ( 1.0_wp - exp_term )
+                particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
+                                       ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
+                                       g / exp_arg ) * ( 1.0_wp - exp_term )
+             ENDIF
+
+          ENDIF
+
+       ENDDO
+    
+    ELSE
+
+       DO  n = 1, number_of_particles
+
+!--       Transport of particles with inertia
+          particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
+          particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
+          particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
+!
+!--       Update of the particle velocity
+          IF ( cloud_droplets )  THEN
+!
+!--          Terminal velocity is computed for vertical direction (Rogers et al., 
+!--          1993, J. Appl. Meteorol.)
+             diameter = particles(n)%radius * 2000.0_wp !diameter in mm
+             IF ( diameter <= d0_rog )  THEN
+                w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
+             ELSE
+                w_s = a_rog - b_rog * EXP( -c_rog * diameter )
+             ENDIF
+
+!
+!--          If selected, add random velocities following Soelch and Kaercher
+!--          (2010, Q. J. R. Meteorol. Soc.)
+             IF ( use_sgs_for_particles )  THEN
+                 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
+                 RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
+                 sigma          = SQRT( e(kp,jp,ip) )
+
+                 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+                 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+                 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+
+                 particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
+                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg1
+                 particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
+                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg2
+                 particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
+                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg3
+
+                 particles(n)%speed_x = u_int(n) + particles(n)%rvar1
+                 particles(n)%speed_y = v_int(n) + particles(n)%rvar2
+                 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
+             ELSE
+                 particles(n)%speed_x = u_int(n)
+                 particles(n)%speed_y = v_int(n)
+                 particles(n)%speed_z = w_int(n) - w_s
+             ENDIF
+
+          ELSE
+
+             IF ( use_sgs_for_particles )  THEN
                 exp_arg  = particle_groups(particles(n)%group)%exp_arg
                 exp_term = EXP( -exp_arg * dt_particle(n) )
@@ -1353 +1470 @@
                 exp_term = particle_groups(particles(n)%group)%exp_term
              ENDIF
-             particles(n)%speed_x = particles(n)%speed_x * exp_term +          &
+             particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
                                     u_int(n) * ( 1.0_wp - exp_term )
-             particles(n)%speed_y = particles(n)%speed_y * exp_term +          &
+             particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
                                     v_int(n) * ( 1.0_wp - exp_term )
-             particles(n)%speed_z = particles(n)%speed_z * exp_term +          &
-                                    ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) *  &
-                                    g / exp_arg ) * ( 1.0_wp - exp_term )
+             particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
+                                    ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
+                                    exp_arg ) * ( 1.0_wp - exp_term )
           ENDIF
 
-       ENDDO
-    
-    ELSE
-
-       DO  n = 1, number_of_particles
-
-!--       Transport of particles with inertia
-          particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
-          particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
-          particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
-!
-!--       Update of the particle velocity
-          IF ( cloud_droplets )  THEN
-
-             exp_arg  = 4.5_wp * dens_ratio(n) * molecular_viscosity /         &
-                        ( particles(n)%radius )**2 *                           &
-                        ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius *  &
-                          SQRT( ( u_int(n) - particles(n)%speed_x )**2 +       &
-                                ( v_int(n) - particles(n)%speed_y )**2 +       &
-                                ( w_int(n) - particles(n)%speed_z )**2 ) /     &
-                                         molecular_viscosity )**0.687_wp       &
-                        )
-
-             exp_term = EXP( -exp_arg * dt_particle(n) )
-          ELSEIF ( use_sgs_for_particles )  THEN
-             exp_arg  = particle_groups(particles(n)%group)%exp_arg
-             exp_term = EXP( -exp_arg * dt_particle(n) )
-          ELSE
-             exp_arg  = particle_groups(particles(n)%group)%exp_arg
-             exp_term = particle_groups(particles(n)%group)%exp_term
-          ENDIF
-          particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
-                                 u_int(n) * ( 1.0_wp - exp_term )
-          particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
-                                 v_int(n) * ( 1.0_wp - exp_term )
-          particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
-                                 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
-                                 exp_arg ) * ( 1.0_wp - exp_term )
        ENDDO
 

palm/trunk/SOURCE/lpm_boundary_conds.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Tails removed. Unused variables removed.
 ! 
 ! Former revisions:
@@ -75 +75 @@
 
     USE control_parameters,                                                    &
-        ONLY:  dz, message_string, particle_maximum_age, simulated_time
+        ONLY:  dz, message_string, particle_maximum_age
 
     USE cpulog,                                                                &
@@ -89 +89 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  deleted_particles, deleted_tails, ibc_par_b, ibc_par_t,         &
-               number_of_particles, particles,                                 &
-               particle_tail_coordinates, particle_type, offset_ocean_nzt_m1,  &
-               tail_mask, use_particle_tails, use_sgs_for_particles
+        ONLY:  deleted_particles, ibc_par_b, ibc_par_t, number_of_particles,   &
+               particles, particle_type, offset_ocean_nzt_m1,                  &
+               use_sgs_for_particles
 
     USE pegrid
@@ -119 +118 @@
     INTEGER(iwp) ::  k5             !<
     INTEGER(iwp) ::  n              !<
-    INTEGER(iwp) ::  nn             !<
     INTEGER(iwp) ::  t_index        !<
     INTEGER(iwp) ::  t_index_number !<
@@ -150 +148 @@
        DO  n = 1, number_of_particles
 
-          nn = particles(n)%tail_id
-
 !
 !--       Stop if particles have moved further than the length of one 
@@ -171 +167 @@
              particles(n)%particle_mask  = .FALSE.
              deleted_particles = deleted_particles + 1
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                tail_mask(nn) = .FALSE.
-                deleted_tails = deleted_tails + 1
-             ENDIF
           ENDIF
 
@@ -183 +175 @@
                 particles(n)%particle_mask  = .FALSE.
                 deleted_particles = deleted_particles + 1
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   tail_mask(nn) = .FALSE.
-                   deleted_tails = deleted_tails + 1
-                ENDIF
              ELSEIF ( ibc_par_t == 2 )  THEN
 !
@@ -196 +184 @@
                    particles(n)%rvar3 = -particles(n)%rvar3
                 ENDIF
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zu(nz) - &
-                                               particle_tail_coordinates(1,3,nn)
-                ENDIF
              ENDIF
           ENDIF
@@ -209 +193 @@
                 particles(n)%particle_mask  = .FALSE.
                 deleted_particles = deleted_particles + 1
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   tail_mask(nn) = .FALSE.
-                   deleted_tails = deleted_tails + 1
-                ENDIF
              ELSEIF ( ibc_par_b == 2 )  THEN
 !
@@ -221 +201 @@
                      particles(n)%rvar3 < 0.0_wp )  THEN
                    particles(n)%rvar3 = -particles(n)%rvar3
-                ENDIF
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zu(nz) - &
-                                               particle_tail_coordinates(1,3,nn)
-                ENDIF
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zw(0) - &
-                                               particle_tail_coordinates(1,3,nn)
                 ENDIF
              ENDIF

palm/trunk/SOURCE/lpm_calc_liquid_water_content.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Unused variables removed.
 ! 
 ! Former revisions:
@@ -85 +85 @@
     INTEGER(iwp) ::  k   !<
     INTEGER(iwp) ::  n   !<
-    INTEGER(iwp) ::  psi !<
 
     CALL cpu_log( log_point_s(45), 'lpm_calc_ql', 'start' )

palm/trunk/SOURCE/lpm_collision_kernels.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-!
+! PALM kernel has been deleted.
+! Bugfix in the calculation of the turbulent enhancement factor of the
+! collection efficiency.
+!
+! Unused variables removed.
 ! 
 ! Former revisions:
@@ -96 +100 @@
 !> This module calculates collision efficiencies either due to pure gravitational
 !> effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
-!> including the effects of (SGS) turbulence (Wang kernel, see Wang and
-!> Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been
+!> including the effects of turbulence (Wang kernel, see Wang and
+!> Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8, and Ayala et al., 2008: 
+!> New J. Phys., 10, 075016). The original code has been
 !> provided by L.-P. Wang but is substantially reformatted and speed optimized
 !> here.
-!>
-!> @attention Physical quantities (like g, densities, etc.) used in this module
-!>            still have to be adjusted to those values used in the main PALM code.
-!>            Also, quantities in CGS-units should be converted to SI-units eventually.
 !------------------------------------------------------------------------------!
  MODULE lpm_collision_kernels_mod
@@ -114 +115 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  collision_kernel, dissipation_classes, particles, radius_classes
+        ONLY:  collision_kernel, dissipation_classes, particles,               &
+               radius_classes
 
     USE pegrid
@@ -123 +125 @@
     PRIVATE
 
-    PUBLIC  ckernel, collision_efficiency_rogers, init_kernels, &
-            rclass_lbound, rclass_ubound, recalculate_kernel
+    PUBLIC  ckernel, init_kernels, rclass_lbound, rclass_ubound,               &
+            recalculate_kernel
 
     REAL(wp) ::  epsilon       !<
-    REAL(wp) ::  eps2          !<
     REAL(wp) ::  rclass_lbound !<
     REAL(wp) ::  rclass_ubound !<
     REAL(wp) ::  urms          !<
-    REAL(wp) ::  urms2         !<
-
-    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  epsclass !<
-    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  radclass !<
-    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  winf     !<
+
+    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  epsclass  !< dissipation rate class
+    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  radclass  !< radius class
+    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  winf      !<
     
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ec       !<
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ecf      !<
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gck      !<
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hkernel  !<
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hwratio  !<
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ec        !<
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ecf       !<
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gck       !<
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hkernel   !<
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hwratio   !<
     
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  ckernel !<
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ckernel !<
 
     SAVE
@@ -149 +149 @@
 !
 !-- Public interfaces
-    INTERFACE collision_efficiency_rogers
-       MODULE PROCEDURE collision_efficiency_rogers
-    END INTERFACE collision_efficiency_rogers
-
     INTERFACE init_kernels
        MODULE PROCEDURE init_kernels
@@ -186 +182 @@
        IF ( collision_kernel(6:9) == 'fast' )  THEN
 
-          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,               &
-                    0:dissipation_classes), epsclass(1:dissipation_classes), &
+          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,                 &
+                    0:dissipation_classes), epsclass(1:dissipation_classes),   &
                     radclass(1:radius_classes) )
 
@@ -195 +191 @@
           rclass_lbound = LOG( 1.0E-6_wp )
           rclass_ubound = LOG( 2.0E-4_wp )
-          radclass(1)   = 1.0E-6_wp
+          radclass(1)   = EXP( rclass_lbound )
           DO  i = 2, radius_classes
              radclass(i) = EXP( rclass_lbound +                                &
@@ -377 +373 @@
 ! Description:
 ! ------------
-!> Calculation of gck
-!> This is from Aayala 2008b, page 37ff.
-!> Necessary input parameters: water density, radii of droplets, air density,
-!> air viscosity, turbulent dissipation rate, taylor microscale reynolds number,
-!> gravitational acceleration  --> to be replaced by PALM parameters
+!> Calculation of effects of turbulence on the geometric collision kernel 
+!> (by including the droplets' average radial relative velocities and their 
+!> radial distribution function) following the analytic model by Aayala et al. 
+!> (2008, New J. Phys.). For details check the second part 2 of the publication,
+!> page 37ff.
+!>
+!> Input parameters, which need to be replaced by PALM parameters:
+!>    water density, air density
 !------------------------------------------------------------------------------!
     SUBROUTINE turbsd
@@ -392 +391 @@
 
        IMPLICIT NONE
-       
-       LOGICAL, SAVE ::  first = .TRUE. !<
 
        INTEGER(iwp) ::  i     !<
@@ -442 +439 @@
        REAL(wp) ::  z         !<
 
-       REAL(wp), DIMENSION(1:radius_classes) ::  st  !<
-       REAL(wp), DIMENSION(1:radius_classes) ::  tau !<
-       
-!
-!--    Initial assignment of constants
-       IF ( first )  THEN
-
-          first = .FALSE.
-
-       ENDIF
-
-       lambda    = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon ) ! in m
+       REAL(wp), DIMENSION(1:radius_classes) ::  st  !< Stokes number
+       REAL(wp), DIMENSION(1:radius_classes) ::  tau !< inertial time scale
+
+       lambda    = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon )
        lambda_re = urms**2 * SQRT( 15.0_wp / epsilon / molecular_viscosity )
-       tl        = urms**2 / epsilon                       ! in s
-       lf        = 0.5_wp * urms**3 / epsilon              ! in m
-       tauk      = SQRT( molecular_viscosity / epsilon )                  ! in s
-       eta       = ( molecular_viscosity**3 / epsilon )**0.25_wp          ! in m
+       tl        = urms**2 / epsilon
+       lf        = 0.5_wp * urms**3 / epsilon
+       tauk      = SQRT( molecular_viscosity / epsilon )
+       eta       = ( molecular_viscosity**3 / epsilon )**0.25_wp
        vk        = eta / tauk
 
        ao = ( 11.0_wp + 7.0_wp * lambda_re ) / ( 205.0_wp + lambda_re )
-       tt = SQRT( 2.0_wp * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk  ! in s
-
-       CALL fallg    ! gives winf in m/s
+       tt = SQRT( 2.0_wp * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk
+
+!
+!--    Get terminal velocity of droplets
+       CALL fallg
 
        DO  i = 1, radius_classes
-          tau(i) = winf(i) / g    ! in s
-          st(i)  = tau(i) / tauk
+          tau(i) = winf(i) / g    ! inertial time scale
+          st(i)  = tau(i) / tauk  ! Stokes number
        ENDDO
 
 !
-!--    Calculate wr (from Aayala 2008b, page 38f)
+!--    Calculate average radial relative velocity at contact (wrfin)
        z   = tt / tl
        be  = SQRT( 2.0_wp ) * lambda / lf
        bbb = SQRT( 1.0_wp - 2.0_wp * be**2 )
        d1  = ( 1.0_wp + bbb ) / ( 2.0_wp * bbb )
-       e1  = lf * ( 1.0_wp + bbb ) * 0.5_wp   ! in m
+       e1  = lf * ( 1.0_wp + bbb ) * 0.5_wp
        d2  = ( 1.0_wp - bbb ) * 0.5_wp / bbb
-       e2  = lf * ( 1.0_wp - bbb ) * 0.5_wp   ! in m
+       e2  = lf * ( 1.0_wp - bbb ) * 0.5_wp
        ccc = SQRT( 1.0_wp - 2.0_wp * z**2 )
        b1  = ( 1.0_wp + ccc ) * 0.5_wp / ccc
-       c1  = tl * ( 1.0_wp + ccc ) * 0.5_wp   ! in s
+       c1  = tl * ( 1.0_wp + ccc ) * 0.5_wp
        b2  = ( 1.0_wp - ccc ) * 0.5_wp / ccc
-       c2  = tl * ( 1.0_wp - ccc ) * 0.5_wp   ! in s
+       c2  = tl * ( 1.0_wp - ccc ) * 0.5_wp
 
        DO  i = 1, radius_classes
 
-          v1 = winf(i)        ! in m/s
-          t1 = tau(i)         ! in s
+          v1 = winf(i)
+          t1 = tau(i)
 
           DO  j = 1, i
              rrp = radclass(i) + radclass(j)
-             v2  = winf(j)                                 ! in m/s
-             t2  = tau(j)                                  ! in s
+             v2  = winf(j)
+             t2  = tau(j)
 
              v1xysq  = b1 * d1 * phi_w(c1,e1,v1,t1) - b1 * d2 * phi_w(c1,e2,v1,t1) &
                      - b2 * d1 * phi_w(c2,e1,v1,t1) + b2 * d2 * phi_w(c2,e2,v1,t1)
-             v1xysq  = v1xysq * urms**2 / t1                ! in m**2/s**2
-             vrms1xy = SQRT( v1xysq )                       ! in m/s
+             v1xysq  = v1xysq * urms**2 / t1
+             vrms1xy = SQRT( v1xysq )
 
              v2xysq  = b1 * d1 * phi_w(c1,e1,v2,t2) - b1 * d2 * phi_w(c1,e2,v2,t2) &
                      - b2 * d1 * phi_w(c2,e1,v2,t2) + b2 * d2 * phi_w(c2,e2,v2,t2)
-             v2xysq  = v2xysq * urms**2 / t2                ! in m**2/s**2
-             vrms2xy = SQRT( v2xysq )                       ! in m/s
+             v2xysq  = v2xysq * urms**2 / t2
+             vrms2xy = SQRT( v2xysq )
 
              IF ( winf(i) >= winf(j) )  THEN
@@ -523 +514 @@
                          b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
              fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
-             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)   ! in m**2/s**2
-             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0_wp * v1v2xy ! in m**2/s**2
+             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)
+             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0_wp * v1v2xy
              IF ( wrtur2xy < 0.0_wp )  wrtur2xy = 0.0_wp
              wrgrav2  = pi / 8.0_wp * ( winf(j) - winf(i) )**2
-             wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) ) ! in m/s
-
-!
-!--          Calculate gr
+             wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) )
+
+!
+!--          Calculate radial distribution function (grfin)
              IF ( st(j) > st(i) )  THEN
                 sst = st(j)
@@ -546 +537 @@
              ao_gr  = ao + ( pi / 8.0_wp) * ( g / vk * tauk )**2
              fao_gr = 20.115_wp * SQRT( ao_gr / lambda_re )
-             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
+             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta
 
              grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5_wp )
              IF ( grfin < 1.0_wp )  grfin = 1.0_wp
 
-             gck(i,j) = 2.0_wp * pi * rrp**2 * wrfin * grfin        ! in cm**3/s
+!
+!--          Calculate general collection kernel (without the consideration of
+!--          collection efficiencies)
+             gck(i,j) = 2.0_wp * pi * rrp**2 * wrfin * grfin
              gck(j,i) = gck(i,j)
 
@@ -559 +553 @@
     END SUBROUTINE turbsd
 
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> phi_w as a function
-!------------------------------------------------------------------------------!
     REAL(wp) FUNCTION phi_w( a, b, vsett, tau0 )
-
+!
+!--    Function used in the Ayala et al. (2008) analytical model for turbulent
+!--    effects on the collision kernel
        IMPLICIT NONE
 
@@ -576 +566 @@
 
        aa1 = 1.0_wp / tau0 + 1.0_wp / a + vsett / b
-       phi_w = 1.0_wp / aa1  - 0.5_wp * vsett / b / aa1**2  ! in s
+       phi_w = 1.0_wp / aa1  - 0.5_wp * vsett / b / aa1**2
 
     END FUNCTION phi_w
 
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> zhi as a function
-!------------------------------------------------------------------------------!
     REAL(wp) FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
-
+!
+!--    Function used in the Ayala et al. (2008) analytical model for turbulent
+!--    effects on the collision kernel
        IMPLICIT NONE
 
@@ -613 +599 @@
              b / aa3**2 + ( 4.0_wp / aa4 - 1.0_wp / aa5**2 - 1.0_wp / aa1**2 ) &
              * vsett2 * 0.5_wp / b /aa6 + ( 2.0_wp * ( b / aa2 - b / aa1 ) -   &
-             vsett1 / aa2**2 + vsett2 / aa1**2 ) * 0.5_wp / b / aa3    ! in s**2
+             vsett1 / aa2**2 + vsett2 / aa1**2 ) * 0.5_wp / b / aa3
 
     END FUNCTION zhi
@@ -621 +607 @@
 ! Description:
 ! ------------
-!> Calculation of terminal velocity winf following Equations 10-138 to 10-145
-!> from (Pruppacher and Klett, 1997)
+!> Parameterization of terminal velocity following Rogers et al. (1993, J. Appl.
+!> Meteorol.)
 !------------------------------------------------------------------------------!
     SUBROUTINE fallg
- 
-       USE cloud_parameters,                                                   &
-           ONLY:  rho_l
-    
-       USE control_parameters,                                                 &
-           ONLY:  g
 
        USE particle_attributes,                                                &
            ONLY:  radius_classes
 
-
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i !<
-       INTEGER(iwp) ::  j !<
-
-       LOGICAL, SAVE ::  first = .TRUE. !<
-
-       REAL(wp), SAVE ::  cunh  !<
-       REAL(wp), SAVE ::  eta   !<
-       REAL(wp), SAVE ::  phy   !<
-       REAL(wp), SAVE ::  py    !<
-       REAL(wp), SAVE ::  rho_a !<
-       REAL(wp), SAVE ::  sigma !<
-       REAL(wp), SAVE ::  stb   !<
-       REAL(wp), SAVE ::  stok  !<
-       REAL(wp), SAVE ::  xlamb !<
-
-       REAL(wp) ::  bond        !<
-       REAL(wp) ::  x           !<
-       REAL(wp) ::  xrey        !<
-       REAL(wp) ::  y           !<
-
-       REAL(wp), DIMENSION(1:7), SAVE  ::  b !<
-       REAL(wp), DIMENSION(1:6), SAVE  ::  c !<
-
-!
-!--    Initial assignment of constants
-       IF ( first )  THEN
-
-          first = .FALSE.
-          b = (/  -0.318657E1_wp,   0.992696E0_wp,  -0.153193E-2_wp, &
-                  -0.987059E-3_wp, -0.578878E-3_wp,  0.855176E-4_wp, &
-                  -0.327815E-5_wp /)
-          c = (/  -0.500015E1_wp,   0.523778E1_wp,  -0.204914E1_wp,   &
-                   0.475294E0_wp,  -0.542819E-1_wp,  0.238449E-2_wp /)
-
-!
-!--       Parameter values for p = 1013,25 hPa and T = 293,15 K
-          eta   = 1.818E-5_wp         ! in kg/(m s)
-          xlamb = 6.6E-8_wp           ! in m
-          rho_a = 1.204_wp            ! in kg/m**3
-          cunh  = 1.26_wp * xlamb     ! in m
-          sigma = 0.07363_wp          ! in kg/s**2
-          stok  = 2.0_wp  * g * ( rho_l - rho_a ) / ( 9.0_wp * eta ) ! in 1/(m s)
-          stb   = 32.0_wp * rho_a * ( rho_l - rho_a) * g / (3.0_wp * eta * eta)
-          phy   = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) )
-          py    = phy**( 1.0_wp / 6.0_wp )
-
-       ENDIF
+       INTEGER(iwp) ::  j                            !<
+
+       REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp    !< parameter
+       REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp   !< parameter
+       REAL(wp), PARAMETER ::  a_rog     = 9.65_wp   !< parameter
+       REAL(wp), PARAMETER ::  b_rog     = 10.43_wp  !< parameter
+       REAL(wp), PARAMETER ::  c_rog     = 0.6_wp    !< parameter
+       REAL(wp), PARAMETER ::  d0_rog    = 0.745_wp  !< seperation diameter
+
+       REAL(wp)            ::  diameter              !< droplet diameter in mm
+
 
        DO  j = 1, radius_classes
 
-          IF ( radclass(j) <= 1.0E-5_wp )  THEN
-
-             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) )
-
-          ELSEIF ( radclass(j) > 1.0E-5_wp  .AND.  radclass(j) <= 5.35E-4_wp )  THEN
-
-             x = LOG( stb * radclass(j)**3 )
-             y = 0.0_wp
-
-             DO  i = 1, 7
-                y = y + b(i) * x**(i-1)
-             ENDDO
-!
-!--          Note: this Eq. is wrong in (Pruppacher and Klett, 1997, p. 418)
-!--          for correct version see (Beard, 1976)
-             xrey = ( 1.0_wp + cunh / radclass(j) ) * EXP( y )
-
-             winf(j) = xrey * eta / ( 2.0_wp * rho_a * radclass(j) )
-
-          ELSEIF ( radclass(j) > 5.35E-4_wp )  THEN
-
-             IF ( radclass(j) > 0.0035_wp )  THEN
-                bond = g * ( rho_l - rho_a ) * 0.0035_wp**2 / sigma
-             ELSE
-               bond = g * ( rho_l - rho_a ) * radclass(j)**2 / sigma
-             ENDIF
-
-             x = LOG( 16.0_wp * bond * py / 3.0_wp )
-             y = 0.0_wp
-
-             DO  i = 1, 6
-                y = y + c(i) * x**(i-1)
-             ENDDO
-
-             xrey = py * EXP( y )
-
-             IF ( radclass(j) > 0.0035_wp )  THEN
-                winf(j) = xrey * eta / ( 2.0_wp * rho_a * 0.0035_wp )
-             ELSE
-                winf(j) = xrey * eta / ( 2.0_wp * rho_a * radclass(j) )
-             ENDIF
-
+          diameter = radclass(j) * 2000.0_wp
+
+          IF ( diameter <= d0_rog )  THEN
+             winf(j) = k_cap_rog * diameter * ( 1.0_wp -                       &
+                                                EXP( -k_low_rog * diameter ) )
+          ELSE
+             winf(j) = a_rog - b_rog * EXP( -c_rog * diameter )
           ENDIF
 
@@ -740 +648 @@
 ! Description:
 ! ------------
-!> Calculation of collision efficiencies for the Hall kernel
+!> Interpolation of collision efficiencies (Hall, 1980, J. Atmos. Sci.)
 !------------------------------------------------------------------------------!
     SUBROUTINE effic
@@ -776 +684 @@
 
          first = .FALSE.
-         r0  = (/   6.0_wp,   8.0_wp,  10.0_wp, 15.0_wp,  20.0_wp,  25.0_wp,  &
-                   30.0_wp,  40.0_wp,  50.0_wp, 60.0_wp,  70.0_wp, 100.0_wp,  &
+         r0  = (/   6.0_wp,   8.0_wp,  10.0_wp, 15.0_wp,  20.0_wp,  25.0_wp,   &
+                   30.0_wp,  40.0_wp,  50.0_wp, 60.0_wp,  70.0_wp, 100.0_wp,   &
                   150.0_wp, 200.0_wp, 300.0_wp /)
 
-         rat = (/ 0.00_wp, 0.05_wp, 0.10_wp, 0.15_wp, 0.20_wp, 0.25_wp,       &
-                  0.30_wp, 0.35_wp, 0.40_wp, 0.45_wp, 0.50_wp, 0.55_wp,       &
-                  0.60_wp, 0.65_wp, 0.70_wp, 0.75_wp, 0.80_wp, 0.85_wp,       &
+         rat = (/ 0.00_wp, 0.05_wp, 0.10_wp, 0.15_wp, 0.20_wp, 0.25_wp,        &
+                  0.30_wp, 0.35_wp, 0.40_wp, 0.45_wp, 0.50_wp, 0.55_wp,        &
+                  0.60_wp, 0.65_wp, 0.70_wp, 0.75_wp, 0.80_wp, 0.85_wp,        &
                   0.90_wp, 0.95_wp, 1.00_wp /)
 
-         ecoll(:,1)  = (/ 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, &
-                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, &
+         ecoll(:,1)  = (/ 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp,    &
+                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp,    &
                           0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp /)
-         ecoll(:,2)  = (/ 0.003_wp, 0.003_wp, 0.003_wp, 0.004_wp, 0.005_wp, &
-                          0.005_wp, 0.005_wp, 0.010_wp, 0.100_wp, 0.050_wp, &
+         ecoll(:,2)  = (/ 0.003_wp, 0.003_wp, 0.003_wp, 0.004_wp, 0.005_wp,    &
+                          0.005_wp, 0.005_wp, 0.010_wp, 0.100_wp, 0.050_wp,    &
                           0.200_wp, 0.500_wp, 0.770_wp, 0.870_wp, 0.970_wp /)
-         ecoll(:,3)  = (/ 0.007_wp, 0.007_wp, 0.007_wp, 0.008_wp, 0.009_wp, &
-                          0.010_wp, 0.010_wp, 0.070_wp, 0.400_wp, 0.430_wp, &
+         ecoll(:,3)  = (/ 0.007_wp, 0.007_wp, 0.007_wp, 0.008_wp, 0.009_wp,    &
+                          0.010_wp, 0.010_wp, 0.070_wp, 0.400_wp, 0.430_wp,    &
                           0.580_wp, 0.790_wp, 0.930_wp, 0.960_wp, 1.000_wp /)
-         ecoll(:,4)  = (/ 0.009_wp, 0.009_wp, 0.009_wp, 0.012_wp, 0.015_wp, &
-                          0.010_wp, 0.020_wp, 0.280_wp, 0.600_wp, 0.640_wp, &
+         ecoll(:,4)  = (/ 0.009_wp, 0.009_wp, 0.009_wp, 0.012_wp, 0.015_wp,    &
+                          0.010_wp, 0.020_wp, 0.280_wp, 0.600_wp, 0.640_wp,    &
                           0.750_wp, 0.910_wp, 0.970_wp, 0.980_wp, 1.000_wp /)
-         ecoll(:,5)  = (/ 0.014_wp, 0.014_wp, 0.014_wp, 0.015_wp, 0.016_wp, &
-                          0.030_wp, 0.060_wp, 0.500_wp, 0.700_wp, 0.770_wp, &
+         ecoll(:,5)  = (/ 0.014_wp, 0.014_wp, 0.014_wp, 0.015_wp, 0.016_wp,    &
+                          0.030_wp, 0.060_wp, 0.500_wp, 0.700_wp, 0.770_wp,    &
                           0.840_wp, 0.950_wp, 0.970_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,6)  = (/ 0.017_wp, 0.017_wp, 0.017_wp, 0.020_wp, 0.022_wp, &
-                          0.060_wp, 0.100_wp, 0.620_wp, 0.780_wp, 0.840_wp, &
+         ecoll(:,6)  = (/ 0.017_wp, 0.017_wp, 0.017_wp, 0.020_wp, 0.022_wp,    &
+                          0.060_wp, 0.100_wp, 0.620_wp, 0.780_wp, 0.840_wp,    &
                           0.880_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,7)  = (/ 0.030_wp, 0.030_wp, 0.024_wp, 0.022_wp, 0.032_wp, &
-                          0.062_wp, 0.200_wp, 0.680_wp, 0.830_wp, 0.870_wp, &
+         ecoll(:,7)  = (/ 0.030_wp, 0.030_wp, 0.024_wp, 0.022_wp, 0.032_wp,    &
+                          0.062_wp, 0.200_wp, 0.680_wp, 0.830_wp, 0.870_wp,    &
                           0.900_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,8)  = (/ 0.025_wp, 0.025_wp, 0.025_wp, 0.036_wp, 0.043_wp, &
-                          0.130_wp, 0.270_wp, 0.740_wp, 0.860_wp, 0.890_wp, &
+         ecoll(:,8)  = (/ 0.025_wp, 0.025_wp, 0.025_wp, 0.036_wp, 0.043_wp,    &
+                          0.130_wp, 0.270_wp, 0.740_wp, 0.860_wp, 0.890_wp,    &
                           0.920_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,9)  = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.040_wp, 0.052_wp, &
-                          0.200_wp, 0.400_wp, 0.780_wp, 0.880_wp, 0.900_wp, &
+         ecoll(:,9)  = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.040_wp, 0.052_wp,    &
+                          0.200_wp, 0.400_wp, 0.780_wp, 0.880_wp, 0.900_wp,    &
                           0.940_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,10) = (/ 0.030_wp, 0.030_wp, 0.030_wp, 0.047_wp, 0.064_wp, &
-                          0.250_wp, 0.500_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
+         ecoll(:,10) = (/ 0.030_wp, 0.030_wp, 0.030_wp, 0.047_wp, 0.064_wp,    &
+                          0.250_wp, 0.500_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
                           0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,11) = (/ 0.040_wp, 0.040_wp, 0.033_wp, 0.037_wp, 0.068_wp, &
-                          0.240_wp, 0.550_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
+         ecoll(:,11) = (/ 0.040_wp, 0.040_wp, 0.033_wp, 0.037_wp, 0.068_wp,    &
+                          0.240_wp, 0.550_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
                           0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,12) = (/ 0.035_wp, 0.035_wp, 0.035_wp, 0.055_wp, 0.079_wp, &
-                          0.290_wp, 0.580_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
+         ecoll(:,12) = (/ 0.035_wp, 0.035_wp, 0.035_wp, 0.055_wp, 0.079_wp,    &
+                          0.290_wp, 0.580_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
                           0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,13) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.062_wp, 0.082_wp, &
-                          0.290_wp, 0.590_wp, 0.780_wp, 0.900_wp, 0.910_wp, &
+         ecoll(:,13) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.062_wp, 0.082_wp,    &
+                          0.290_wp, 0.590_wp, 0.780_wp, 0.900_wp, 0.910_wp,    &
                           0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,14) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.060_wp, 0.080_wp, &
-                          0.290_wp, 0.580_wp, 0.770_wp, 0.890_wp, 0.910_wp, &
+         ecoll(:,14) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.060_wp, 0.080_wp,    &
+                          0.290_wp, 0.580_wp, 0.770_wp, 0.890_wp, 0.910_wp,    &
                           0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,15) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.041_wp, 0.075_wp, &
-                          0.250_wp, 0.540_wp, 0.760_wp, 0.880_wp, 0.920_wp, &
+         ecoll(:,15) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.041_wp, 0.075_wp,    &
+                          0.250_wp, 0.540_wp, 0.760_wp, 0.880_wp, 0.920_wp,    &
                           0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,16) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.052_wp, 0.067_wp, &
-                          0.250_wp, 0.510_wp, 0.770_wp, 0.880_wp, 0.930_wp, &
+         ecoll(:,16) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.052_wp, 0.067_wp,    &
+                          0.250_wp, 0.510_wp, 0.770_wp, 0.880_wp, 0.930_wp,    &
                           0.970_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,17) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.047_wp, 0.057_wp, &
-                          0.250_wp, 0.490_wp, 0.770_wp, 0.890_wp, 0.950_wp, &
+         ecoll(:,17) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.047_wp, 0.057_wp,    &
+                          0.250_wp, 0.490_wp, 0.770_wp, 0.890_wp, 0.950_wp,    &
                           1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
-         ecoll(:,18) = (/ 0.036_wp, 0.036_wp, 0.036_wp, 0.042_wp, 0.048_wp, &
-                          0.230_wp, 0.470_wp, 0.780_wp, 0.920_wp, 1.000_wp, &
+         ecoll(:,18) = (/ 0.036_wp, 0.036_wp, 0.036_wp, 0.042_wp, 0.048_wp,    &
+                          0.230_wp, 0.470_wp, 0.780_wp, 0.920_wp, 1.000_wp,    &
                           1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp /)
-         ecoll(:,19) = (/ 0.040_wp, 0.040_wp, 0.035_wp, 0.033_wp, 0.040_wp, &
-                          0.112_wp, 0.450_wp, 0.790_wp, 1.010_wp, 1.030_wp, &
+         ecoll(:,19) = (/ 0.040_wp, 0.040_wp, 0.035_wp, 0.033_wp, 0.040_wp,    &
+                          0.112_wp, 0.450_wp, 0.790_wp, 1.010_wp, 1.030_wp,    &
                           1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp /)
-         ecoll(:,20) = (/ 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, &
-                          0.119_wp, 0.470_wp, 0.950_wp, 1.300_wp, 1.700_wp, &
+         ecoll(:,20) = (/ 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp,    &
+                          0.119_wp, 0.470_wp, 0.950_wp, 1.300_wp, 1.700_wp,    &
                           2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp /)
-         ecoll(:,21) = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, &
-                          0.125_wp, 0.520_wp, 1.400_wp, 2.300_wp, 3.000_wp, &
+         ecoll(:,21) = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp,    &
+                          0.125_wp, 0.520_wp, 1.400_wp, 2.300_wp, 3.000_wp,    &
                           4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp /)
        ENDIF
@@ -852 +760 @@
 !
 !--    Calculate the radius class index of particles with respect to array r
-!--    Radius has to be in µm
+!--    Radius has to be in microns
        ALLOCATE( ira(1:radius_classes) )
        DO  j = 1, radius_classes
@@ -867 +775 @@
 !
 !--    Two-dimensional linear interpolation of the collision efficiency.
-!--    Radius has to be in µm
+!--    Radius has to be in microns
        DO  j = 1, radius_classes
           DO  i = 1, j
@@ -878 +786 @@
              IF ( ir < 16 )  THEN
                 IF ( ir >= 2 )  THEN
-                   pp = ( ( MAX(radclass(j),radclass(i)) * 1.0E06_wp ) - r0(ir-1) ) / &
-                        ( r0(ir) - r0(ir-1) )
+                   pp = ( ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp ) -     &
+                          r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
                    qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
                    ec(j,i) = ( 1.0_wp - pp ) * ( 1.0_wp - qq )                 &
@@ -911 +819 @@
 ! Description:
 ! ------------
-!> Calculation of enhancement factor for collision efficencies due to turbulence
+!> Interpolation of turbulent enhancement factor for collision efficencies
+!> following Wang and Grabowski (2009, Atmos. Sci. Let.)
 !------------------------------------------------------------------------------!
     SUBROUTINE turb_enhance_eff
@@ -940 +849 @@
 
        REAL(wp), DIMENSION(1:11), SAVE ::  rat           !<
-       
        REAL(wp), DIMENSION(1:7), SAVE  ::  r0            !<
        
@@ -958 +866 @@
                    0.7_wp, 0.8_wp, 0.9_wp, 1.0_wp /)
 !
-!--       for 100 cm**2/s**3
+!--       Tabulated turbulent enhancement factor at 100 cm**2/s**3
           ecoll_100(:,1)  = (/  1.74_wp,   1.74_wp,   1.773_wp, 1.49_wp,  &
                                 1.207_wp,  1.207_wp,  1.0_wp /)
@@ -979 +887 @@
           ecoll_100(:,10) = (/  1.570_wp,  1.570_wp,  1.244_wp, 1.166_wp, &
                                 1.088_wp,  1.088_wp,  1.0_wp /)
-          ecoll_100(:,11) = (/ 20.3_wp,   20.3_wp,   14.6_wp,  8.61_wp,  &
+          ecoll_100(:,11) = (/ 20.3_wp,   20.3_wp,   14.6_wp,   8.61_wp,  &
                                 2.60_wp,   2.60_wp,   1.0_wp /)
 !
-!--       for 400 cm**2/s**3
+!--       Tabulated turbulent enhancement factor at 400 cm**2/s**3
           ecoll_400(:,1)  = (/  4.976_wp,  4.976_wp,  3.593_wp,  2.519_wp, &
                                 1.445_wp,  1.445_wp,  1.0_wp /)
@@ -1010 +918 @@
 !
 !--    Calculate the radius class index of particles with respect to array r0
-!--    Radius has to be in µm
+!--    The droplet radius has to be given in microns.
        ALLOCATE( ira(1:radius_classes) )
 
@@ -1025 +933 @@
 
 !
-!--    Two-dimensional linear interpolation of the collision efficiencies
-!--    Radius has to be in µm
+!--    Two-dimensional linear interpolation of the turbulent enhancement factor.
+!--    The droplet radius has to be given in microns.
        DO  j =  1, radius_classes
           DO  i = 1, j
@@ -1040 +948 @@
              ENDDO
 
-             y1 = 0.0001_wp      ! for 0 m**2/s**3
+             y1 = 1.0_wp  ! turbulent enhancement factor at 0 m**2/s**3
 
              IF ( ir < 8 )  THEN
                 IF ( ir >= 2 )  THEN
-                   pp = ( MAX(radclass(j),radclass(i))*1.0E6_wp - r0(ir-1) ) / &
-                        ( r0(ir) - r0(ir-1) )
+                   pp = ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp -  &
+                          r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
                    qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
                    y2 = ( 1.0_wp - pp ) * ( 1.0_wp - qq ) * ecoll_100(ir-1,iq-1) + &
@@ -1066 +974 @@
              ENDIF
 !
-!--          Linear interpolation of dissipation rate in m**2/s**3
+!--          Linear interpolation of turbulent enhancement factor
              IF ( epsilon <= 0.01_wp )  THEN
                 ecf(j,i) = ( epsilon - 0.01_wp ) / ( 0.0_wp  - 0.01_wp ) * y1 &
@@ -1087 +995 @@
     END SUBROUTINE turb_enhance_eff
 
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Collision efficiencies from table 8.2 in Rogers and Yau (1989, 3rd edition).
-!> Values are calculated from table by bilinear interpolation.
-!------------------------------------------------------------------------------!
- 
-    SUBROUTINE collision_efficiency_rogers( mean_r, r, e)
-
-       IMPLICIT NONE
-
-       INTEGER(iwp)  ::  i !<
-       INTEGER(iwp)  ::  j !<
-       INTEGER(iwp)  ::  k !<
-
-       LOGICAL, SAVE ::  first = .TRUE. !<
-
-       REAL(wp)      ::  aa      !<
-       REAL(wp)      ::  bb      !<
-       REAL(wp)      ::  cc      !<
-       REAL(wp)      ::  dd      !<
-       REAL(wp)      ::  dx      !<
-       REAL(wp)      ::  dy      !<
-       REAL(wp)      ::  e       !<
-       REAL(wp)      ::  gg      !<
-       REAL(wp)      ::  mean_r  !<
-       REAL(wp)      ::  mean_rm !<
-       REAL(wp)      ::  r       !<
-       REAL(wp)      ::  rm      !<
-       REAL(wp)      ::  x       !<
-       REAL(wp)      ::  y       !<
- 
-       REAL(wp), DIMENSION(1:9), SAVE      ::  collected_r = 0.0_wp !<
-       
-       REAL(wp), DIMENSION(1:19), SAVE     ::  collector_r = 0.0_wp !<
-       
-       REAL(wp), DIMENSION(1:9,1:19), SAVE ::  ef = 0.0_wp          !<
-
-       mean_rm = mean_r * 1.0E06_wp
-       rm      = r      * 1.0E06_wp
-
-       IF ( first )  THEN
-
-          collected_r = (/    2.0_wp,    3.0_wp,    4.0_wp,    6.0_wp,    8.0_wp, &
-                             10.0_wp,   15.0_wp,   20.0_wp,   25.0_wp /)
-          collector_r = (/   10.0_wp,   20.0_wp,   30.0_wp,   40.0_wp,   50.0_wp, &
-                             60.0_wp,   80.0_wp,  100.0_wp,  150.0_wp,  200.0_wp, &
-                            300.0_wp,  400.0_wp,  500.0_wp,  600.0_wp, 1000.0_wp, &
-                           1400.0_wp, 1800.0_wp, 2400.0_wp, 3000.0_wp /)
-
-          ef(:,1)  = (/ 0.017_wp, 0.027_wp, 0.037_wp, 0.052_wp, 0.052_wp,      &
-                        0.052_wp, 0.052_wp, 0.0_wp,   0.0_wp /)
-          ef(:,2)  = (/ 0.001_wp, 0.016_wp, 0.027_wp, 0.060_wp, 0.12_wp,       &
-                        0.17_wp,  0.17_wp,  0.17_wp,  0.0_wp /)
-          ef(:,3)  = (/ 0.001_wp, 0.001_wp, 0.02_wp,  0.13_wp,  0.28_wp,       &
-                        0.37_wp,  0.54_wp,  0.55_wp,  0.47_wp/)
-          ef(:,4)  = (/ 0.001_wp, 0.001_wp, 0.02_wp,  0.23_wp,  0.4_wp,        &
-                        0.55_wp,  0.7_wp,   0.75_wp,  0.75_wp/)
-          ef(:,5)  = (/ 0.01_wp,  0.01_wp,  0.03_wp,  0.3_wp,   0.4_wp,        &
-                        0.58_wp,  0.73_wp,  0.75_wp,  0.79_wp/)
-          ef(:,6)  = (/ 0.01_wp,  0.01_wp,  0.13_wp,  0.38_wp,  0.57_wp,       &
-                        0.68_wp,  0.80_wp,  0.86_wp,  0.91_wp/)
-          ef(:,7)  = (/ 0.01_wp,  0.085_wp, 0.23_wp,  0.52_wp,  0.68_wp,       &
-                        0.76_wp,  0.86_wp,  0.92_wp,  0.95_wp/)
-          ef(:,8)  = (/ 0.01_wp,  0.14_wp,  0.32_wp,  0.60_wp,  0.73_wp,       &
-                        0.81_wp,  0.90_wp,  0.94_wp,  0.96_wp/)
-          ef(:,9)  = (/ 0.025_wp, 0.25_wp,  0.43_wp,  0.66_wp,  0.78_wp,       &
-                        0.83_wp,  0.92_wp,  0.95_wp,  0.96_wp/)
-          ef(:,10) = (/ 0.039_wp, 0.3_wp,   0.46_wp,  0.69_wp,  0.81_wp,       &
-                        0.87_wp,  0.93_wp,  0.95_wp,  0.96_wp/)
-          ef(:,11) = (/ 0.095_wp, 0.33_wp,  0.51_wp,  0.72_wp,  0.82_wp,       &
-                        0.87_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
-          ef(:,12) = (/ 0.098_wp, 0.36_wp,  0.51_wp,  0.73_wp,  0.83_wp,       &
-                        0.88_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
-          ef(:,13) = (/ 0.1_wp,   0.36_wp,  0.52_wp,  0.74_wp,  0.83_wp,       &
-                        0.88_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
-          ef(:,14) = (/ 0.17_wp,  0.4_wp,   0.54_wp,  0.72_wp,  0.83_wp,       &
-                        0.88_wp,  0.94_wp,  0.98_wp,  1.0_wp /)
-          ef(:,15) = (/ 0.15_wp,  0.37_wp,  0.52_wp,  0.74_wp,  0.82_wp,       &
-                        0.88_wp,  0.94_wp,  0.98_wp,  1.0_wp /)
-          ef(:,16) = (/ 0.11_wp,  0.34_wp,  0.49_wp,  0.71_wp,  0.83_wp,       &
-                        0.88_wp,  0.94_wp,  0.95_wp,  1.0_wp /)
-          ef(:,17) = (/ 0.08_wp,  0.29_wp,  0.45_wp,  0.68_wp,  0.8_wp,        &
-                        0.86_wp,  0.96_wp,  0.94_wp,  1.0_wp /)
-          ef(:,18) = (/ 0.04_wp,  0.22_wp,  0.39_wp,  0.62_wp,  0.75_wp,       &
-                        0.83_wp,  0.92_wp,  0.96_wp,  1.0_wp /)
-          ef(:,19) = (/ 0.02_wp,  0.16_wp,  0.33_wp,  0.55_wp,  0.71_wp,       &
-                        0.81_wp,  0.90_wp,  0.94_wp,  1.0_wp /)
-
-       ENDIF
-
-       DO  k = 1, 8
-          IF ( collected_r(k) <= mean_rm )  i = k
-       ENDDO
-
-       DO  k = 1, 18
-          IF ( collector_r(k) <= rm )  j = k
-       ENDDO
-
-       IF ( rm < 10.0_wp )  THEN
-          e = 0.0_wp
-       ELSEIF ( mean_rm < 2.0_wp )  THEN
-          e = 0.001_wp
-       ELSEIF ( mean_rm >= 25.0_wp )  THEN
-          IF( j <= 2 )  e = 0.0_wp
-          IF( j == 3 )  e = 0.47_wp
-          IF( j == 4 )  e = 0.8_wp
-          IF( j == 5 )  e = 0.9_wp
-          IF( j >=6  )  e = 1.0_wp
-       ELSEIF ( rm >= 3000.0_wp )  THEN
-          IF( i == 1 )  e = 0.02_wp
-          IF( i == 2 )  e = 0.16_wp
-          IF( i == 3 )  e = 0.33_wp
-          IF( i == 4 )  e = 0.55_wp
-          IF( i == 5 )  e = 0.71_wp
-          IF( i == 6 )  e = 0.81_wp
-          IF( i == 7 )  e = 0.90_wp
-          IF( i >= 8 )  e = 0.94_wp
-       ELSE
-          x  = mean_rm - collected_r(i)
-          y  = rm - collector_r(j)
-          dx = collected_r(i+1) - collected_r(i)
-          dy = collector_r(j+1) - collector_r(j)
-          aa = x**2 + y**2
-          bb = ( dx - x )**2 + y**2
-          cc = x**2 + ( dy - y )**2
-          dd = ( dx - x )**2 + ( dy - y )**2
-          gg = aa + bb + cc + dd
-
-          e = ( (gg-aa)*ef(i,j) + (gg-bb)*ef(i+1,j) + (gg-cc)*ef(i,j+1) + &
-                (gg-dd)*ef(i+1,j+1) ) / (3.0_wp * gg)
-       ENDIF
-
-    END SUBROUTINE collision_efficiency_rogers
-
  END MODULE lpm_collision_kernels_mod

palm/trunk/SOURCE/lpm_data_output_particles.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Tails removed. Unused variables removed.
 !
 ! Former revisions:
@@ -58 +58 @@
 
     USE control_parameters,                                                    &
-        ONLY:  prt_time_count, simulated_time
+        ONLY:  simulated_time
 
     USE cpulog,                                                                &
@@ -68 +68 @@
     USE kinds
 
-    USE netcdf_interface,                                                      &
-        ONLY:  nc_stat
-
     USE particle_attributes,                                                   &
-        ONLY:  grid_particles, maximum_number_of_particles,                    &
-               maximum_number_of_tailpoints, maximum_number_of_tails,          &
-               number_of_particles, number_of_tails, particles,                &
-               particle_tail_coordinates, prt_count
+        ONLY:  grid_particles, number_of_particles,  particles, prt_count
 
     IMPLICIT NONE
@@ -103 +97 @@
        ENDDO
     ENDDO
-!
-!-- particle tails currently not available
-!    WRITE ( 85 )  maximum_number_of_tailpoints, maximum_number_of_tails, &
-!                  number_of_tails
-!    IF ( maximum_number_of_tails > 0 )  THEN
-!       WRITE ( 85 )  particle_tail_coordinates, prt_time_count
-!    ENDIF
 
     CALL close_file( 85 )

palm/trunk/SOURCE/lpm_droplet_collision.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Integration of a new collision algortithm based on Shima et al. (2009) and
+! Soelch and Kaercher (2010) called all_or_nothing. The previous implemented 
+! collision algorithm is called average_impact. Moreover, both algorithms are 
+! now positive definit due to their construction, i.e., no negative weighting 
+! factors should occur.
 ! 
 ! Former revisions:
@@ -63 +67 @@
 !> Calculates change in droplet radius by collision. Droplet collision is
 !> calculated for each grid box seperately. Collision is parameterized by
-!> using collision kernels. Three different kernels are available:
-!> PALM kernel: Kernel is approximated using a method from Rogers and
-!>              Yau (1989, A Short Course in Cloud Physics, Pergamon Press).
-!>              All droplets smaller than the treated one are represented by
-!>              one droplet with mean features. Collision efficiencies are taken
-!>              from the respective table in Rogers and Yau.
+!> using collision kernels. Two different kernels are available:
 !> Hall kernel: Kernel from Hall (1980, J. Atmos. Sci., 2486-2507), which
 !>              considers collision due to pure gravitational effects.
@@ -85 +84 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  diss, ql, ql_v, ql_vp, u, v, w, zu, zw
+        ONLY:  diss, ql_v, ql_vp
 
     USE cloud_parameters,                                                      &
-        ONLY:  effective_coll_efficiency
+        ONLY:  rho_l
 
     USE constants,                                                             &
@@ -94 +93 @@
 
     USE control_parameters,                                                    &
-        ONLY:  dt_3d, message_string, u_gtrans, v_gtrans, dz
+        ONLY:  dt_3d, message_string, dz
 
     USE cpulog,                                                                &
@@ -100 +99 @@
 
     USE grid_variables,                                                        &
-        ONLY:  ddx, dx, ddy, dy
-
-    USE indices,                                                               &
-        ONLY:  nxl, nxr, nyn, nys, nzb, nzt 
+        ONLY:  dx, dy
 
     USE kinds
 
     USE lpm_collision_kernels_mod,                                             &
-        ONLY:  ckernel, collision_efficiency_rogers, recalculate_kernel
+        ONLY:  ckernel, recalculate_kernel
 
     USE particle_attributes,                                                   &
-        ONLY:  deleted_particles, dissipation_classes, hall_kernel,            &
-               palm_kernel, particles, particle_type,                          &
-               prt_count, use_kernel_tables, wang_kernel
+        ONLY:  all_or_nothing, average_impact, dissipation_classes,            &
+               hall_kernel, iran_part, number_of_particles, particles,         &
+               particle_type, prt_count, use_kernel_tables, wang_kernel
+
+    USE random_function_mod,                                                   &
+        ONLY:  random_function
 
     USE pegrid
@@ -121 +120 @@
     INTEGER(iwp) ::  eclass   !<
     INTEGER(iwp) ::  i        !<
-    INTEGER(iwp) ::  ii       !<
-    INTEGER(iwp) ::  inc      !<
-    INTEGER(iwp) ::  is       !<
     INTEGER(iwp) ::  j        !<
-    INTEGER(iwp) ::  jj       !<
-    INTEGER(iwp) ::  js       !<
     INTEGER(iwp) ::  k        !<
-    INTEGER(iwp) ::  kk       !<
     INTEGER(iwp) ::  n        !<
-    INTEGER(iwp) ::  pse      !<
-    INTEGER(iwp) ::  psi      !<
+    INTEGER(iwp) ::  m        !<
     INTEGER(iwp) ::  rclass_l !<
     INTEGER(iwp) ::  rclass_s !<
 
-    INTEGER(iwp), DIMENSION(prt_count(k,j,i)) ::  rclass_v !<
-
-    LOGICAL, SAVE ::  first_flag = .TRUE. !<
-
-    TYPE(particle_type) :: tmp_particle   !<
-
-    REAL(wp) ::  aa       !<
-    REAL(wp) ::  auxn     !< temporary variables
-    REAL(wp) ::  auxs     !< temporary variables
-    REAL(wp) ::  bb       !<
-    REAL(wp) ::  cc       !<
-    REAL(wp) ::  dd       !<
-    REAL(wp) ::  ddV      !<
-    REAL(wp) ::  delta_r  !<
-    REAL(wp) ::  delta_v  !<
-    REAL(wp) ::  epsilon  !<
-    REAL(wp) ::  gg       !<
-    REAL(wp) ::  mean_r   !<
-    REAL(wp) ::  ql_int   !<
-    REAL(wp) ::  ql_int_l !<
-    REAL(wp) ::  ql_int_u !<
-    REAL(wp) ::  r3       !<
-    REAL(wp) ::  sl_r3    !<
-    REAL(wp) ::  sl_r4    !<
-    REAL(wp) ::  sum1     !<
-    REAL(wp) ::  sum2     !<
-    REAL(wp) ::  sum3     !<
-    REAL(wp) ::  u_int    !<
-    REAL(wp) ::  u_int_l  !<
-    REAL(wp) ::  u_int_u  !<
-    REAL(wp) ::  v_int    !<
-    REAL(wp) ::  v_int_l  !<
-    REAL(wp) ::  v_int_u  !<
-    REAL(wp) ::  w_int    !<
-    REAL(wp) ::  w_int_l  !<
-    REAL(wp) ::  w_int_u  !<
-    REAL(wp) ::  x        !<
-    REAL(wp) ::  y        !<
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad    !<
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight !<
-
-    REAL, DIMENSION(prt_count(k,j,i))    :: ck
-    REAL, DIMENSION(prt_count(k,j,i))    :: r3v
-    REAL, DIMENSION(prt_count(k,j,i))    :: sum1v
-    REAL, DIMENSION(prt_count(k,j,i))    :: sum2v
+    REAL(wp) ::  collection_probability  !< probability for collection
+    REAL(wp) ::  ddV                     !< inverse grid box volume
+    REAL(wp) ::  epsilon                 !< dissipation rate
+    REAL(wp) ::  factor_volume_to_mass   !< 4.0 / 3.0 * pi * rho_l
+    REAL(wp) ::  xm                      !< mean mass of droplet m
+    REAL(wp) ::  xn                      !< mean mass of droplet n
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight  !< weighting factor
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mass    !< total mass of super droplet
 
     CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
 
-!
-!-- Collision requires at least two particles in the box
-    IF ( prt_count(k,j,i) > 1 )  THEN
-!
-!--    First, sort particles within the gridbox by their size,
-!--    using Shell's method (see Numerical Recipes)
-!--    NOTE: In case of using particle tails, the re-sorting of
-!--    ----  tails would have to be included here!
-       IF ( .NOT. ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
-       use_kernel_tables ) )  THEN
-          psi = 0
-          inc = 1
-          DO WHILE ( inc <= prt_count(k,j,i) )
-             inc = 3 * inc + 1
-          ENDDO
-
-          DO WHILE ( inc > 1 )
-             inc = inc / 3
-             DO  is = inc+1, prt_count(k,j,i)
-                tmp_particle = particles(psi+is)
-                js = is
-                DO WHILE ( particles(psi+js-inc)%radius >             &
-                tmp_particle%radius )
-                   particles(psi+js) = particles(psi+js-inc)
-                   js = js - inc
-                   IF ( js <= inc )  EXIT
-                ENDDO
-                particles(psi+js) = tmp_particle
-             ENDDO
-          ENDDO
-       ENDIF
-
-       psi = 1
-       pse = prt_count(k,j,i)
+    number_of_particles   = prt_count(k,j,i)
+    factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l 
+    ddV                   = 1 / ( dx * dy * dz )
+!
+!-- Collision requires at least one super droplet inside the box
+    IF ( number_of_particles > 0 )  THEN
 
 !
 !--    Now apply the different kernels
-       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
-            use_kernel_tables )  THEN
-!
-!--       Fast method with pre-calculated efficiencies for
+       IF ( use_kernel_tables )  THEN
+!
+!--       Fast method with pre-calculated collection kernels for
 !--       discrete radius- and dissipation-classes.
 !--
@@ -244 +170 @@
 !--       Droplet collision are calculated using collision-coalescence
 !--       formulation proposed by Wang (see PALM documentation)
-!--       Since new radii after collision are defined by radii of all
-!--       droplets before collision, temporary fields for new radii and
-!--       weighting factors are needed
-          ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
-
-          rad    = 0.0_wp
-          weight = 0.0_wp
-
-          sum1v(1:prt_count(k,j,i)) = 0.0_wp
-          sum2v(1:prt_count(k,j,i)) = 0.0_wp
-
-          DO  n = 1, prt_count(k,j,i)
-
-             rclass_l = particles(n)%class
-!
-!--          Mass added due to collisions with smaller droplets
-             DO  is = n+1, prt_count(k,j,i)
-                rclass_s = particles(is)%class
-                auxs = ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor
-                auxn = ckernel(rclass_s,rclass_l,eclass) * particles(n)%weight_factor
-                IF ( particles(is)%radius <  particles(n)%radius )  THEN
-                   sum1v(n) =  sum1v(n)  + particles(is)%radius**3 * auxs
-                   sum2v(is) = sum2v(is) + auxn
-                ELSE
-                   sum2v(n)  = sum2v(n)  + auxs
-                   sum1v(is) = sum1v(is) + particles(n)%radius**3 * auxn
-                ENDIF
+!--       Temporary fields for total mass of super-droplet and weighting factors 
+!--       are allocated.
+          ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
+
+          mass(1:number_of_particles)   = particles(1:number_of_particles)%weight_factor * &
+                                          particles(1:number_of_particles)%radius**3     * &
+                                          factor_volume_to_mass
+          weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
+
+          IF ( average_impact )  THEN  ! select collision algorithm
+
+             DO  n = 1, number_of_particles
+
+                rclass_l = particles(n)%class
+                xn       = mass(n) / weight(n)
+
+                DO  m = n, number_of_particles
+
+                   rclass_s = particles(m)%class
+                   xm       = mass(m) / weight(m)
+
+                   IF ( xm .LT. xn )  THEN
+                      
+!
+!--                   Particle n collects smaller particle m
+                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
+                                               weight(n) * ddV * dt_3d
+
+                      mass(n)   = mass(n)   + mass(m)   * collection_probability
+                      weight(m) = weight(m) - weight(m) * collection_probability
+                      mass(m)   = mass(m)   - mass(m)   * collection_probability
+                   ELSEIF ( xm .GT. xn )  THEN 
+!
+!--                   Particle m collects smaller particle n
+                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
+                                               weight(m) * ddV * dt_3d
+
+                      mass(m)   = mass(m)   + mass(n)   * collection_probability
+                      weight(n) = weight(n) - weight(n) * collection_probability
+                      mass(n)   = mass(n)   - mass(n)   * collection_probability
+                   ELSE
+!
+!--                   Same-size collections. If n = m, weight is reduced by the
+!--                   number of possible same-size collections; the total mass 
+!--                   is not changed during same-size collection.
+!--                   Same-size collections of different 
+!--                   particles ( n /= m ) are treated as same-size collections 
+!--                   of ONE partilce with weight = weight(n) + weight(m) and 
+!--                   mass = mass(n) + mass(m).
+!--                   Accordingly, each particle loses the same number of 
+!--                   droplets to the other particle, but this has no effect on 
+!--                   total mass mass, since the exchanged droplets have the 
+!--                   same radius. 
+
+!--                   Note: For m = n this equation is an approximation only 
+!--                   valid for weight >> 1 (which is usually the case). The
+!--                   approximation is weight(n)-1 = weight(n).
+                      weight(n) = weight(n) - 0.5_wp * weight(n) *                &
+                                              ckernel(rclass_l,rclass_s,eclass) * &
+                                              weight(m) * ddV * dt_3d
+                      IF ( n .NE. m )  THEN 
+                         weight(m) = weight(m) - 0.5_wp * weight(m) *                &
+                                                 ckernel(rclass_l,rclass_s,eclass) * &
+                                                 weight(n) * ddV * dt_3d
+                      ENDIF
+                   ENDIF
+
+                ENDDO
+
+                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
+
              ENDDO
-          ENDDO
-          rclass_v = particles(1:prt_count(k,j,i))%class
-          DO  n = 1, prt_count(k,j,i)
-            ck(n)       = ckernel(rclass_v(n),rclass_v(n),eclass)
-          ENDDO
-          r3v      = particles(1:prt_count(k,j,i))%radius**3
-          DO  n = 1, prt_count(k,j,i)
-             sum3 = 0.0_wp
-             ddV = ddx * ddy / dz
-!
-!--          Change of the current weighting factor
-             sum3 = 1 - dt_3d * ddV *                                 &
-                        ck(n) *                                       &
-                        ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
-                    dt_3d * ddV * sum2v(n)
-             weight(n) = particles(n)%weight_factor * sum3
-!
-!--          Change of the current droplet radius
-             rad(n) = ( (r3v(n) + dt_3d * ddV * (sum1v(n) - sum2v(n) * r3v(n)) )/&
-                              sum3 )**0.33333333333333_wp
-
-             ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n) &
-                                         * rad(n)**3
-
-          ENDDO
+
+          ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
+
+             DO  n = 1, number_of_particles
+
+                rclass_l = particles(n)%class
+                xn       = mass(n) / weight(n) ! mean mass of droplet n
+
+                DO  m = n, number_of_particles
+
+                   rclass_s = particles(m)%class
+                   xm = mass(m) / weight(m) ! mean mass of droplet m
+
+                   IF ( weight(n) .LT. weight(m) )  THEN
+!
+!--                   Particle n collects weight(n) droplets of particle m  
+                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
+                                               weight(m) * ddV * dt_3d
+
+                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
+                         mass(n)   = mass(n)   + weight(n) * xm
+                         weight(m) = weight(m) - weight(n)
+                         mass(m)   = mass(m)   - weight(n) * xm
+                      ENDIF
+
+                   ELSEIF ( weight(m) .LT. weight(n) )  THEN 
+!
+!--                   Particle m collects weight(m) droplets of particle n  
+                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
+                                               weight(n) * ddV * dt_3d
+
+                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
+                         mass(m)   = mass(m)   + weight(m) * xn
+                         weight(n) = weight(n) - weight(m)
+                         mass(n)   = mass(n)   - weight(m) * xn
+                      ENDIF
+                   ELSE
+!
+!--                   Collisions of particles of the same weighting factor. 
+!--                   Particle n collects 1/2 weight(n) droplets of particle m,
+!--                   particle m collects 1/2 weight(m) droplets of particle n.
+!--                   The total mass mass changes accordingly.
+!--                   If n = m, the first half of the droplets coalesces with the 
+!--                   second half of the droplets; mass is unchanged because
+!--                   xm = xn for n = m.
+
+!--                   Note: For m = n this equation is an approximation only 
+!--                   valid for weight >> 1 (which is usually the case). The
+!--                   approximation is weight(n)-1 = weight(n).
+                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
+                                               weight(n) * ddV * dt_3d
+
+                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
+                         mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
+                         mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
+                         weight(n) = weight(n) - 0.5_wp * weight(m)
+                         weight(m) = weight(n)
+                      ENDIF
+                   ENDIF
+
+                ENDDO
+
+                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
+
+             ENDDO
+
+          ENDIF
+
+
+
+
           IF ( ANY(weight < 0.0_wp) )  THEN
                 WRITE( message_string, * ) 'negative weighting'
@@ -303 +319 @@
           ENDIF
 
-          particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
-          particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
-
-          DEALLOCATE(rad, weight)
-
-       ELSEIF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
-                .NOT. use_kernel_tables )  THEN
-!
-!--       Collision efficiencies are calculated for every new
+          particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
+                                                      ( weight(1:number_of_particles) &
+                                                        * factor_volume_to_mass       &
+                                                      )                               &
+                                                    )**0.33333333333333_wp
+
+          particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
+
+          DEALLOCATE(weight, mass)
+
+       ELSEIF ( .NOT. use_kernel_tables )  THEN
+!
+!--       Collection kernels are calculated for every new
 !--       grid box. First, allocate memory for kernel table.
 !--       Third dimension is 1, because table is re-calculated for
 !--       every new dissipation value.
-          ALLOCATE( ckernel(1:prt_count(k,j,i),1:prt_count(k,j,i),1:1) )
-!
-!--       Now calculate collision efficiencies for this box
+          ALLOCATE( ckernel(1:number_of_particles,1:number_of_particles,1:1) )
+!
+!--       Now calculate collection kernel for this box. Note that
+!--       the kernel is based on the previous time step
           CALL recalculate_kernel( i, j, k )
-
 !
 !--       Droplet collision are calculated using collision-coalescence
 !--       formulation proposed by Wang (see PALM documentation)
-!--       Since new radii after collision are defined by radii of all
-!--       droplets before collision, temporary fields for new radii and
-!--       weighting factors are needed
-          ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
-
-          rad = 0.0_wp
-          weight = 0.0_wp
-
-          DO  n = psi, pse, 1
-
-             sum1 = 0.0_wp
-             sum2 = 0.0_wp
-             sum3 = 0.0_wp
-!
-!--          Mass added due to collisions with smaller droplets
-             DO  is = psi, n-1
-                sum1 = sum1 + ( particles(is)%radius**3 *            &
-                                ckernel(n,is,1) *  &
-                                particles(is)%weight_factor )
+!--       Temporary fields for total mass of super-droplet and weighting factors 
+!--       are allocated.
+          ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
+
+          mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
+                                        particles(1:number_of_particles)%radius**3     * &
+                                        factor_volume_to_mass
+
+          weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
+
+          IF ( average_impact )  THEN  ! select collision algorithm
+
+             DO  n = 1, number_of_particles
+
+                xn = mass(n) / weight(n) ! mean mass of droplet n
+
+                DO  m = n, number_of_particles
+
+                   xm = mass(m) / weight(m) !mean mass of droplet m
+
+                   IF ( xm .LT. xn )  THEN
+!
+!--                   Particle n collects smaller particle m 
+                      collection_probability = ckernel(n,m,1) * weight(n) *    &
+                                               ddV * dt_3d
+
+                      mass(n)   = mass(n)   + mass(m)   * collection_probability
+                      weight(m) = weight(m) - weight(m) * collection_probability
+                      mass(m)   = mass(m)   - mass(m)   * collection_probability
+                   ELSEIF ( xm .GT. xn )  THEN 
+!
+!--                   Particle m collects smaller particle n
+                      collection_probability = ckernel(n,m,1) * weight(m) *    &
+                                               ddV * dt_3d
+
+                      mass(m)   = mass(m)   + mass(n)   * collection_probability
+                      weight(n) = weight(n) - weight(n) * collection_probability
+                      mass(n)   = mass(n)   - mass(n)   * collection_probability
+                   ELSE
+!
+!--                   Same-size collections. If n = m, weight is reduced by the
+!--                   number of possible same-size collections; the total mass 
+!--                   mass is not changed during same-size collection. 
+!--                   Same-size collections of different 
+!--                   particles ( n /= m ) are treated as same-size collections 
+!--                   of ONE partilce with weight = weight(n) + weight(m) and 
+!--                   mass = mass(n) + mass(m).
+!--                   Accordingly, each particle loses the same number of 
+!--                   droplets to the other particle, but this has no effect on 
+!--                   total mass mass, since the exchanged droplets have the 
+!--                   same radius. 
+!--
+!--                   Note: For m = n this equation is an approximation only 
+!--                   valid for weight >> 1 (which is usually the case). The
+!--                   approximation is weight(n)-1 = weight(n).
+                      weight(n) = weight(n) - 0.5_wp * weight(n) *             &
+                                              ckernel(n,m,1) * weight(m) *     &
+                                              ddV * dt_3d
+                      IF ( n .NE. m )  THEN
+                         weight(m) = weight(m) - 0.5_wp * weight(m) *          &
+                                                 ckernel(n,m,1) * weight(n) *  &
+                                                 ddV * dt_3d
+                      ENDIF
+                   ENDIF
+
+
+                ENDDO
+
+                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
+
              ENDDO
-!
-!--          Rate of collisions with larger droplets
-             DO  is = n+1, pse
-                sum2 = sum2 + ( ckernel(n,is,1) *  &
-                                particles(is)%weight_factor )
+
+          ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
+
+             DO  n = 1, number_of_particles
+
+                xn = mass(n) / weight(n) ! mean mass of droplet n
+
+                DO  m = n, number_of_particles
+
+                   xm = mass(m) / weight(m) !mean mass of droplet m
+
+                   IF ( weight(n) .LT. weight(m) )  THEN
+!
+!--                   Particle n collects smaller particle m 
+                      collection_probability = ckernel(n,m,1) * weight(m) *    &
+                                               ddV * dt_3d
+
+                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
+                         mass(n) = mass(n) + weight(n) * xm  
+                         weight(m)    = weight(m)    - weight(n)
+                         mass(m) = mass(m) - weight(n) * xm
+                      ENDIF
+
+                   ELSEIF ( weight(m) .LT. weight(n) )  THEN
+!
+!--                   Particle m collects smaller particle n
+                      collection_probability = ckernel(n,m,1) * weight(n) *    &
+                                               ddV * dt_3d
+
+                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
+                         mass(m) = mass(m) + weight(m) * xn
+                         weight(n)    = weight(n)    - weight(m)
+                         mass(n) = mass(n) - weight(m) * xn
+                      ENDIF
+                   ELSE
+!
+!--                   Collisions of particles of the same weighting factor. 
+!--                   Particle n collects 1/2 weight(n) droplets of particle m,
+!--                   particle m collects 1/2 weight(m) droplets of particle n.
+!--                   The total mass mass changes accordingly.
+!--                   If n = m, the first half of the droplets coalesces with the 
+!--                   second half of the droplets; mass is unchanged because
+!--                   xm = xn for n = m.
+!--
+!--                   Note: For m = n this equation is an approximation only 
+!--                   valid for weight >> 1 (which is usually the case). The
+!--                   approximation is weight(n)-1 = weight(n).
+                      collection_probability = ckernel(n,m,1) * weight(n) *    &
+                                               ddV * dt_3d
+
+                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
+                         mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
+                         mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
+                         weight(n) = weight(n) - 0.5_wp * weight(m)
+                         weight(m) = weight(n)
+                      ENDIF
+                   ENDIF
+
+
+                ENDDO
+
+                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
+
              ENDDO
 
-             r3 = particles(n)%radius**3
-             ddV = ddx * ddy / dz 
-             is = 1
-!
-!--                   Change of the current weighting factor
-             sum3 = 1 - dt_3d * ddV *                                 &
-                        ckernel(n,n,1) *           &
-                        ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
-                    dt_3d * ddV * sum2
-             weight(n-is+1) = particles(n)%weight_factor * sum3
-!
-!--                   Change of the current droplet radius
-             rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
-                              sum3 )**0.33333333333333_wp
-
-             IF ( weight(n-is+1) < 0.0_wp )  THEN
-                WRITE( message_string, * ) 'negative weighting',      &
-                                           'factor: ', weight(n-is+1)
-                CALL message( 'lpm_droplet_collision', 'PA0037',      &
+          ENDIF
+
+          IF ( ANY(weight < 0.0_wp) )  THEN
+                WRITE( message_string, * ) 'negative weighting'
+                CALL message( 'lpm_droplet_collision', 'PA0028',      &
                                2, 2, -1, 6, 1 )
-             ENDIF
-
-             ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
-                                         * rad(n-is+1)**3
-
-          ENDDO
-
-          particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
-          particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
-
-          DEALLOCATE( rad, weight, ckernel )
-
-       ELSEIF ( palm_kernel )  THEN
-!
-!--       PALM collision kernel
-!
-!--       Calculate the mean radius of all those particles which
-!--       are of smaller size than the current particle and
-!--       use this radius for calculating the collision efficiency
-          DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
-
-             sl_r3 = 0.0_wp
-             sl_r4 = 0.0_wp
-
-             DO is = n-1, psi, -1
-                IF ( particles(is)%radius < particles(n)%radius )  &
-                THEN
-                   sl_r3 = sl_r3 + particles(is)%weight_factor     &
-                                   * particles(is)%radius**3
-                   sl_r4 = sl_r4 + particles(is)%weight_factor     &
-                                   * particles(is)%radius**4
-                ENDIF
-             ENDDO
-
-             IF ( ( sl_r3 ) > 0.0_wp )  THEN
-                mean_r = ( sl_r4 ) / ( sl_r3 )
-
-                CALL collision_efficiency_rogers( mean_r,             &
-                                           particles(n)%radius,       &
-                                           effective_coll_efficiency )
-
-             ELSE
-                effective_coll_efficiency = 0.0_wp
-             ENDIF
-
-             IF ( effective_coll_efficiency > 1.0_wp  .OR.            &
-                  effective_coll_efficiency < 0.0_wp )                &
-             THEN
-                WRITE( message_string, * )  'collision_efficien' , &
-                          'cy out of range:' ,effective_coll_efficiency
-                CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
-                              2, -1, 6, 1 )
-             ENDIF
-
-!
-!--          Interpolation of liquid water content
-             ii = particles(n)%x * ddx
-             jj = particles(n)%y * ddy
-             kk = ( particles(n)%z + 0.5_wp * dz ) / dz
-
-             x  = particles(n)%x - ii * dx
-             y  = particles(n)%y - jj * dy
-             aa = x**2          + y**2
-             bb = ( dx - x )**2 + y**2
-             cc = x**2          + ( dy - y )**2
-             dd = ( dx - x )**2 + ( dy - y )**2
-             gg = aa + bb + cc + dd
-
-             ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
-                                                  ql(kk,jj,ii+1)   &
-                        + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
-                                                  ql(kk,jj+1,ii+1) &
-                        ) / ( 3.0_wp * gg )
-
-             ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
-                                                ql(kk+1,jj,ii+1)   &
-                        + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
-                                                ql(kk+1,jj+1,ii+1) &
-                        ) / ( 3.0_wp * gg )
-
-             ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
-                                 ( ql_int_u - ql_int_l )
-
-!
-!--          Interpolate u velocity-component
-             ii = ( particles(n)%x + 0.5_wp * dx ) * ddx
-             jj =   particles(n)%y * ddy
-             kk = ( particles(n)%z + 0.5_wp * dz ) / dz ! only if equidistant
-
-             IF ( ( particles(n)%z - zu(kk) ) > ( 0.5_wp * dz ) ) kk = kk+1
-
-             x  = particles(n)%x + ( 0.5_wp - ii ) * dx
-             y  = particles(n)%y - jj * dy
-             aa = x**2          + y**2
-             bb = ( dx - x )**2 + y**2
-             cc = x**2          + ( dy - y )**2
-             dd = ( dx - x )**2 + ( dy - y )**2
-             gg = aa + bb + cc + dd
-
-             u_int_l = ( (gg-aa) * u(kk,jj,ii)   + (gg-bb) *       &
-                                                  u(kk,jj,ii+1)    &
-                       + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
-                                                  u(kk,jj+1,ii+1)  &
-                       ) / ( 3.0_wp * gg ) - u_gtrans
-             IF ( kk+1 == nzt+1 )  THEN
-                u_int = u_int_l
-             ELSE
-                u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
-                                                 u(kk+1,jj,ii+1)   &
-                          + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
-                                                 u(kk+1,jj+1,ii+1) &
-                          ) / ( 3.0_wp * gg ) - u_gtrans
-                u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
-                                  * ( u_int_u - u_int_l )
-             ENDIF
-
-!
-!--          Same procedure for interpolation of the v velocity-component 
-!--          (adopt index k from u velocity-component)
-             ii =   particles(n)%x * ddx
-             jj = ( particles(n)%y + 0.5_wp * dy ) * ddy
-
-             x  = particles(n)%x - ii * dx
-             y  = particles(n)%y + ( 0.5_wp - jj ) * dy
-             aa = x**2          + y**2
-             bb = ( dx - x )**2 + y**2
-             cc = x**2          + ( dy - y )**2
-             dd = ( dx - x )**2 + ( dy - y )**2
-             gg = aa + bb + cc + dd
-
-             v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *   &
-                                                   v(kk,jj,ii+1)   &
-                       + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
-                                                   v(kk,jj+1,ii+1) &
-                       ) / ( 3.0_wp * gg ) - v_gtrans
-             IF ( kk+1 == nzt+1 )  THEN
-                v_int = v_int_l
-             ELSE
-                v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
-                                                   v(kk+1,jj,ii+1) &
-                          + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
-                                                 v(kk+1,jj+1,ii+1) &
-                          ) / ( 3.0_wp * gg ) - v_gtrans
-                v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
-                                  * ( v_int_u - v_int_l )
-             ENDIF
-
-!
-!--          Same procedure for interpolation of the w velocity-component 
-!--          (adopt index i from v velocity-component)
-             jj = particles(n)%y * ddy
-             kk = particles(n)%z / dz
-
-             x  = particles(n)%x - ii * dx
-             y  = particles(n)%y - jj * dy
-             aa = x**2          + y**2
-             bb = ( dx - x )**2 + y**2
-             cc = x**2          + ( dy - y )**2
-             dd = ( dx - x )**2 + ( dy - y )**2
-             gg = aa + bb + cc + dd
-
-             w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *   &
-                                                     w(kk,jj,ii+1) &
-                       + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
-                                                   w(kk,jj+1,ii+1) &
-                       ) / ( 3.0_wp * gg )
-             IF ( kk+1 == nzt+1 )  THEN
-                w_int = w_int_l
-             ELSE
-                w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
-                                                   w(kk+1,jj,ii+1) &
-                          + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
-                                                 w(kk+1,jj+1,ii+1) &
-                          ) / ( 3.0_wp * gg )
-                w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
-                                  * ( w_int_u - w_int_l )
-             ENDIF
-
-!
-!--          Change in radius due to collision
-             delta_r = effective_coll_efficiency / 3.0_wp          &
-                       * pi * sl_r3 * ddx * ddy / dz               &
-                       * SQRT( ( u_int - particles(n)%speed_x )**2 &
-                             + ( v_int - particles(n)%speed_y )**2 &
-                             + ( w_int - particles(n)%speed_z )**2 &
-                             ) * dt_3d
-!
-!--          Change in volume due to collision
-             delta_v = particles(n)%weight_factor                  &
-                       * ( ( particles(n)%radius + delta_r )**3    &
-                           - particles(n)%radius**3 )
-
-!
-!--          Check if collected particles provide enough LWC for
-!--          volume change of collector particle
-             IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
-
-                delta_r = ( ( sl_r3/particles(n)%weight_factor )               &
-                            + particles(n)%radius**3 )**( 1.0_wp / 3.0_wp )    &
-                            - particles(n)%radius
-
-                DO  is = n-1, psi, -1
-                   IF ( particles(is)%radius < particles(n)%radius )  THEN
-                      particles(is)%weight_factor = 0.0_wp
-                      particles(is)%particle_mask = .FALSE.
-                      deleted_particles = deleted_particles + 1
-                   ENDIF
-                ENDDO
-
-             ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
-
-                DO  is = n-1, psi, -1
-                   IF ( particles(is)%radius < particles(n)%radius &
-                        .AND.  sl_r3 > 0.0_wp )  THEN
-                      particles(is)%weight_factor =                &
-                                   ( ( particles(is)%weight_factor &
-                                   * ( particles(is)%radius**3 ) ) &
-                                   - ( delta_v                     &
-                                   * particles(is)%weight_factor   &
-                                   * ( particles(is)%radius**3 )   &
-                                   / sl_r3 ) )                     &
-                                   / ( particles(is)%radius**3 )
-
-                      IF ( particles(is)%weight_factor < 0.0_wp )  THEN
-                         WRITE( message_string, * ) 'negative ',   &
-                                         'weighting factor: ',     &
-                                         particles(is)%weight_factor
-                         CALL message( 'lpm_droplet_collision', &
-                                       'PA0039',                &
-                                       2, 2, -1, 6, 1 )
-                      ENDIF
-                   ENDIF
-                ENDDO
-
-             ENDIF
-
-             particles(n)%radius = particles(n)%radius + delta_r
-             ql_vp(k,j,i) = ql_vp(k,j,i) +               &
-                            particles(n)%weight_factor * &
-                            ( particles(n)%radius**3 )
-          ENDDO
-
-          ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
-                                        * particles(psi)%radius**3
-
-       ENDIF  ! collision kernel
-
-    ELSE IF ( prt_count(k,j,i) == 1 )  THEN
-
-       psi = 1
-
-!
-!--    Calculate change of weighting factor due to self collision
-       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
-            use_kernel_tables )  THEN
-
-          IF ( wang_kernel )  THEN
-             eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
-                           dissipation_classes ) + 1
-             epsilon = diss(k,j,i)
-          ELSE
-             epsilon = 0.0_wp
           ENDIF
-          IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
-             eclass = 0   ! Hall kernel is used
-          ELSE
-             eclass = MIN( dissipation_classes, eclass )
-          ENDIF
-
-          ddV = ddx * ddy / dz 
-          rclass_l = particles(psi)%class
-          sum3 = 1 - dt_3d * ddV *                                 &
-                     ( ckernel(rclass_l,rclass_l,eclass) *         &
-                       ( particles(psi)%weight_factor-1 ) * 0.5_wp )
-
-          particles(psi)%radius = ( particles(psi)%radius**3 / &
-                                    sum3 )**0.33333333333333_wp
-          particles(psi)%weight_factor = particles(psi)%weight_factor &
-                                         * sum3
-
-       ELSE IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
-                  .NOT. use_kernel_tables )  THEN
-!
-!--       Collision efficiencies are calculated for every new
-!--       grid box. First, allocate memory for kernel table.
-!--       Third dimension is 1, because table is re-calculated for
-!--       every new dissipation value.
-          ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
-!
-!--       Now calculate collision efficiencies for this box
-          CALL recalculate_kernel( i, j, k )
-
-          ddV = ddx * ddy / dz 
-          sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
-                       ( particles(psi)%weight_factor - 1 ) * 0.5_wp ) 
-
-          particles(psi)%radius = ( particles(psi)%radius**3 / &
-                                    sum3 )**0.33333333333333_wp
-          particles(psi)%weight_factor = particles(psi)%weight_factor &
-                                         * sum3
-
-          DEALLOCATE( ckernel )
-       ENDIF 
-
-      ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
-                       particles(psi)%radius**3
+
+          particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
+                                                      ( weight(1:number_of_particles) &
+                                                        * factor_volume_to_mass       &
+                                                      )                               &
+                                                    )**0.33333333333333_wp
+
+          particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
+
+          DEALLOCATE( weight, mass, ckernel )
+
+       ENDIF
+ 
     ENDIF
-
-!
-!-- Check if condensation of LWC was conserved during collision process
+    
+
+!
+!-- Check if LWC is conserved during collision process
     IF ( ql_v(k,j,i) /= 0.0_wp )  THEN
-       IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.             &
+       IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.                      &
             ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp )  THEN
-          WRITE( message_string, * ) 'LWC is not conserved during',&
-                                     ' collision! ',               &
-                                     'LWC after condensation: ',   &
-                                     ql_v(k,j,i),                  &
-                                     ' LWC after collision: ',     &
-                                     ql_vp(k,j,i)
-          CALL message( 'lpm_droplet_collision', 'PA0040',         &
-                        2, 2, -1, 6, 1 )
+          WRITE( message_string, * ) ' LWC is not conserved during',           &
+                                     ' collision! ',                           &
+                                     ' LWC after condensation: ', ql_v(k,j,i), &
+                                     ' LWC after collision: ', ql_vp(k,j,i)
+          CALL message( 'lpm_droplet_collision', 'PA0040', 2, 2, -1, 6, 1 )
        ENDIF
     ENDIF

palm/trunk/SOURCE/lpm_droplet_condensation.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Unused variables removed.
 ! 
 ! Former revisions:
@@ -87 +87 @@
 
     USE cloud_parameters,                                                      &
-        ONLY:  bfactor, curvature_solution_effects, diff_coeff_l,              &
-               eps_ros, l_d_rv, l_v, rho_l,  r_v, thermal_conductivity_l
+        ONLY:  bfactor, curvature_solution_effects, eps_ros, l_d_rv, l_v,      &
+               rho_l,  r_v
 
     USE constants,                                                             &
@@ -94 +94 @@
 
     USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, dt_3d, dz, message_string,                    &
-               molecular_viscosity, rho_surface
+        ONLY:  dt_3d, dz, message_string, molecular_viscosity, rho_surface
+
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point_s
 
     USE grid_variables,                                                        &
-        ONLY:  dx, ddx, dy, ddy
+        ONLY:  dx, dy
 
     USE lpm_collision_kernels_mod,                                             &
@@ -109 +109 @@
     USE particle_attributes,                                                   &
         ONLY:  block_offset, grid_particles, hall_kernel, number_of_particles, &
-               offset_ocean_nzt, offset_ocean_nzt_m1, particles,               &
-               radius_classes, use_kernel_tables, wang_kernel
+               particles, radius_classes, use_kernel_tables, wang_kernel
 
 

palm/trunk/SOURCE/lpm_exchange_horiz.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Tails removed. Unused variables removed.
 !
 ! Former revisions:
@@ -67 +67 @@
 ! Description:
 ! ------------
-!> Exchange of particles (and tails) between the subdomains.
+! Exchange of particles between the subdomains.
 !------------------------------------------------------------------------------!
  MODULE lpm_exchange_horiz_mod
@@ -93 +93 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  alloc_factor, deleted_particles, deleted_tails, grid_particles, &
-               ibc_par_lr, ibc_par_ns, maximum_number_of_tails,                &
-               maximum_number_of_tailpoints, min_nr_particle,                  &
-               mpi_particle_type, number_of_tails, number_of_particles,        &
-               offset_ocean_nzt, particles,                                    &
-               particle_tail_coordinates, particle_type, prt_count,            &
-               tail_mask, trlp_count_sum,                                      &
+        ONLY:  alloc_factor, deleted_particles, grid_particles,                &
+               ibc_par_lr, ibc_par_ns, min_nr_particle,                        &
+               mpi_particle_type, number_of_particles,                         &
+               offset_ocean_nzt, offset_ocean_nzt_m1, particles,               &
+               particle_type, prt_count, trlp_count_sum,                       &
                trlp_count_recv_sum, trnp_count_sum, trnp_count_recv_sum,       &
                trrp_count_sum, trrp_count_recv_sum, trsp_count_sum,            &
-               trsp_count_recv_sum, use_particle_tails, zero_particle
+               trsp_count_recv_sum, zero_particle
 
     USE pegrid
@@ -156 +154 @@
     INTEGER(iwp) ::  j                                       !<
     INTEGER(iwp) ::  jp                                      !<
-    INTEGER(iwp) ::  k                                       !<
     INTEGER(iwp) ::  kp                                      !<
     INTEGER(iwp) ::  n                                       !<
-    INTEGER(iwp) ::  nn                                      !<
-    INTEGER(iwp) ::  tlength                                 !<
     INTEGER(iwp) ::  trlp_count                              !<
     INTEGER(iwp) ::  trlp_count_recv                         !<
@@ -178 +173 @@
     INTEGER(iwp) ::  trspt_count_recv                        !<
 
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trlpt        !<
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trnpt        !<
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trrpt        !<
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trspt        !<
-
     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !<
     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !<
@@ -196 +186 @@
 #if defined( __parallel )
 
+!
+!-- Exchange between subdomains.
+!-- As soon as one particle has moved beyond the boundary of the domain, it
+!-- is included in the relevant transfer arrays and marked for subsequent
+!-- deletion on this PE.
+!-- First sweep for crossings in x direction. Find out first the number of
+!-- particles to be transferred and allocate temporary arrays needed to store
+!-- them.
+!-- For a one-dimensional decomposition along y, no transfer is necessary,
+!-- because the particle remains on the PE, but the particle coordinate has to
+!-- be adjusted.
     trlp_count  = 0
     trlpt_count = 0
@@ -226 +227 @@
                       IF ( i < nxl )  THEN
                          trlp_count = trlp_count + 1
-                         IF ( particles(n)%tail_id /= 0 )  trlpt_count = trlpt_count + 1
                       ELSEIF ( i > nxr )  THEN
                          trrp_count = trrp_count + 1
-                         IF ( particles(n)%tail_id /= 0 )  trrpt_count = trrpt_count + 1
                       ENDIF
                    ENDIF
@@ -247 +246 @@
        trlp = zero_particle
        trrp = zero_particle
-
-       IF ( use_particle_tails )  THEN
-          ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
-                    trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
-          tlength = maximum_number_of_tailpoints * 5
-       ENDIF
 
        trlp_count  = 0
@@ -269 +262 @@
              particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
              DO  n = 1, number_of_particles
-
-                nn = particles(n)%tail_id
 !
 !--             Only those particles that have not been marked as 'deleted' may
@@ -292 +283 @@
                                particles(n)%origin_x = ( nx + 1 ) * dx + &
                                particles(n)%origin_x
-                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                                  i  = particles(n)%tailpoints
-                                  particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
-                                  + particle_tail_coordinates(1:i,1,nn)
-                               ENDIF
                             ELSE
                                trlp_count = trlp_count + 1
@@ -312 +298 @@
                                ENDIF
 
-                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                                  trlpt_count = trlpt_count + 1
-                                  trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
-                                  trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
-                                  trlpt(:,1,trlpt_count)
-                                  tail_mask(nn) = .FALSE.
-                                  deleted_tails = deleted_tails + 1
-                               ENDIF
                             ENDIF
 
@@ -327 +305 @@
                             particles(n)%particle_mask = .FALSE.
                             deleted_particles = deleted_particles + 1
-                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                               tail_mask(nn) = .FALSE.
-                               deleted_tails = deleted_tails + 1
-                            ENDIF
 
                          ELSEIF ( ibc_par_lr == 2 )  THEN
@@ -348 +322 @@
                          deleted_particles = deleted_particles + 1
 
-                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                            trlpt_count = trlpt_count + 1
-                            trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
-                            tail_mask(nn) = .FALSE.
-                            deleted_tails = deleted_tails + 1
-                         ENDIF
                       ENDIF
 
@@ -367 +335 @@
                                particles(n)%origin_x = particles(n)%origin_x - &
                                ( nx + 1 ) * dx
-                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                                  i = particles(n)%tailpoints
-                                  particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
-                                  + particle_tail_coordinates(1:i,1,nn)
-                               ENDIF
                             ELSE
                                trrp_count = trrp_count + 1
@@ -381 +344 @@
                                deleted_particles = deleted_particles + 1
 
-                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                                  trrpt_count = trrpt_count + 1
-                                  trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
-                                  trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
-                                  ( nx + 1 ) * dx
-                                  tail_mask(nn) = .FALSE.
-                                  deleted_tails = deleted_tails + 1
-                               ENDIF
                             ENDIF
 
@@ -396 +351 @@
                             particles(n)%particle_mask = .FALSE.
                             deleted_particles = deleted_particles + 1
-                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                               tail_mask(nn) = .FALSE.
-                               deleted_tails = deleted_tails + 1
-                            ENDIF
 
                          ELSEIF ( ibc_par_lr == 2 )  THEN
@@ -417 +368 @@
                          deleted_particles = deleted_particles + 1
 
-                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                            trrpt_count = trrpt_count + 1
-                            trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
-                            tail_mask(nn) = .FALSE.
-                            deleted_tails = deleted_tails + 1
-                         ENDIF
                       ENDIF
 
@@ -452 +397 @@
        DEALLOCATE(rvrp)
 
-       IF ( use_particle_tails )  THEN
-
-          CALL MPI_SENDRECV( trlpt_count,      1, MPI_INTEGER, pleft,  0, &
-                             trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
-                             comm2d, status, ierr )
-
-          IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
-          THEN
-             IF ( netcdf_data_format < 3 )  THEN
-                message_string = 'maximum_number_of_tails ' //   &
-                                 'needs to be increased ' //     &
-                                 '&but this is not allowed wi'// &
-                                 'th netcdf_data_format < 3'
-                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
-             ELSE
-                CALL lpm_extend_tail_array( trrpt_count_recv )
-             ENDIF
-          ENDIF
-
-          CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL,      &
-                             pleft, 1,                                         &
-                             particle_tail_coordinates(1,1,number_of_tails+1), &
-                             trrpt_count_recv*tlength, MPI_REAL, pright, 1,    &
-                             comm2d, status, ierr )
-!
-!--       Update the tail ids for the transferred particles
-          nn = number_of_tails
-          DO  n = number_of_particles+1, number_of_particles+trrp_count_recv
-             IF ( particles(n)%tail_id /= 0 )  THEN
-                nn = nn + 1
-                particles(n)%tail_id = nn
-             ENDIF
-          ENDDO
-
-       ENDIF
-
 !
 !--    Send right boundary, receive left boundary
@@ -503 +412 @@
        IF ( trlp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvlp(1:trlp_count_recv))
 
-       DEALLOCATE(rvlp)
-
-       IF ( use_particle_tails )  THEN
-
-          CALL MPI_SENDRECV( trrpt_count,      1, MPI_INTEGER, pright, 0, &
-                             trlpt_count_recv, 1, MPI_INTEGER, pleft,  0, &
-                             comm2d, status, ierr )
-
-          IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
-          THEN
-             IF ( netcdf_data_format < 3 )  THEN
-                message_string = 'maximum_number_of_tails ' //   &
-                                 'needs to be increased ' //     &
-                                 '&but this is not allowed wi'// &
-                                 'th netcdf_data_format < 3'
-                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 ) 
-             ELSE
-                CALL lpm_extend_tail_array( trlpt_count_recv )
-             ENDIF
-          ENDIF
-
-          CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL,      &
-                             pright, 1,                                        &
-                             particle_tail_coordinates(1,1,number_of_tails+1), &
-                             trlpt_count_recv*tlength, MPI_REAL, pleft, 1,     &
-                             comm2d, status, ierr )
-!
-!--       Update the tail ids for the transferred particles
-          nn = number_of_tails
-          DO  n = number_of_particles+1, number_of_particles+trlp_count_recv
-             IF ( particles(n)%tail_id /= 0 )  THEN
-                nn = nn + 1
-                particles(n)%tail_id = nn
-             ENDIF
-          ENDDO
-
-       ENDIF
-
-!       number_of_particles = number_of_particles + trlp_count_recv
-!       number_of_tails     = number_of_tails     + trlpt_count_recv
-
-       IF ( use_particle_tails )  THEN
-          DEALLOCATE( trlpt, trrpt )
-       ENDIF
+       DEALLOCATE( rvlp )
        DEALLOCATE( trlp, trrp )
 
@@ -589 +455 @@
                       IF ( j < nys )  THEN
                          trsp_count = trsp_count + 1
-                         IF ( particles(n)%tail_id /= 0 )  trspt_count = trspt_count + 1
                       ELSEIF ( j > nyn )  THEN
                          trnp_count = trnp_count + 1
-                         IF ( particles(n)%tail_id /= 0 )  trnpt_count = trnpt_count + 1
                       ENDIF
                    ENDIF
@@ -609 +473 @@
        trsp = zero_particle
        trnp = zero_particle
-
-       IF ( use_particle_tails )  THEN
-          ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
-                    trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
-          tlength = maximum_number_of_tailpoints * 5
-       ENDIF
 
        trsp_count  = nr_move_south
@@ -633 +491 @@
              particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
              DO  n = 1, number_of_particles
-
-                nn = particles(n)%tail_id
 !
 !--             Only those particles that have not been marked as 'deleted' may
@@ -655 +511 @@
                                particles(n)%origin_y = ( ny + 1 ) * dy + &
                                particles(n)%origin_y
-                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                                  i = particles(n)%tailpoints
-                                  particle_tail_coordinates(1:i,2,nn) =        &
-                                     ( ny+1 ) * dy + particle_tail_coordinates(1:i,2,nn)
-                               ENDIF
                             ELSE
                                trsp_count = trsp_count + 1
@@ -677 +528 @@
                                ENDIF
 
-                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                                  trspt_count = trspt_count + 1
-                                  trspt(:,:,trspt_count) = &
-                                  particle_tail_coordinates(:,:,nn)
-                                  trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
-                                  trspt(:,2,trspt_count)
-                                  tail_mask(nn) = .FALSE.
-                                  deleted_tails = deleted_tails + 1
-                               ENDIF
                             ENDIF
 
@@ -693 +535 @@
                             particles(n)%particle_mask = .FALSE.
                             deleted_particles = deleted_particles + 1
-                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                               tail_mask(nn) = .FALSE.
-                               deleted_tails = deleted_tails + 1
-                            ENDIF
 
                          ELSEIF ( ibc_par_ns == 2 )  THEN
@@ -714 +552 @@
                          deleted_particles = deleted_particles + 1
 
-                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                            trspt_count = trspt_count + 1
-                            trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
-                            tail_mask(nn) = .FALSE.
-                            deleted_tails = deleted_tails + 1
-                         ENDIF
                       ENDIF
 
@@ -733 +565 @@
                                particles(n)%origin_y =                         &
                                   particles(n)%origin_y - ( ny + 1 ) * dy
-                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                                  i = particles(n)%tailpoints
-                                  particle_tail_coordinates(1:i,2,nn) =        &
-                                     - (ny+1) * dy + particle_tail_coordinates(1:i,2,nn)
-                               ENDIF
                             ELSE
                                trnp_count = trnp_count + 1
@@ -747 +574 @@
                                particles(n)%particle_mask = .FALSE.
                                deleted_particles = deleted_particles + 1
-                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                                  trnpt_count = trnpt_count + 1
-                                  trnpt(:,:,trnpt_count) =                     &
-                                     particle_tail_coordinates(:,:,nn)
-                                  trnpt(:,2,trnpt_count) =                     &
-                                     trnpt(:,2,trnpt_count) - ( ny + 1 ) * dy
-                                  tail_mask(nn) = .FALSE.
-                                  deleted_tails = deleted_tails + 1
-                               ENDIF
                             ENDIF
 
@@ -763 +581 @@
                             particles(n)%particle_mask = .FALSE.
                             deleted_particles = deleted_particles + 1
-                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                               tail_mask(nn) = .FALSE.
-                               deleted_tails = deleted_tails + 1
-                            ENDIF
 
                          ELSEIF ( ibc_par_ns == 2 )  THEN
@@ -784 +598 @@
                          deleted_particles = deleted_particles + 1
 
-                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                            trnpt_count = trnpt_count + 1
-                            trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
-                            tail_mask(nn) = .FALSE.
-                            deleted_tails = deleted_tails + 1
-                         ENDIF
                       ENDIF
 
@@ -819 +627 @@
        DEALLOCATE(rvnp)
 
-       IF ( use_particle_tails )  THEN
-
-          CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, psouth, 0, &
-                             trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
-                             comm2d, status, ierr )
-
-          IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
-          THEN
-             IF ( netcdf_data_format < 3 )  THEN
-                message_string = 'maximum_number_of_tails ' //    &
-                                 'needs to be increased ' //      &
-                                 '&but this is not allowed wi' // &
-                                 'th netcdf_data_format < 3'
-                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 ) 
-             ELSE
-                CALL lpm_extend_tail_array( trnpt_count_recv )
-             ENDIF
-          ENDIF
-
-          CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL,      &
-                             psouth, 1,                                        &
-                             particle_tail_coordinates(1,1,number_of_tails+1), &
-                             trnpt_count_recv*tlength, MPI_REAL, pnorth, 1,    &
-                             comm2d, status, ierr )
-
-!
-!--       Update the tail ids for the transferred particles
-          nn = number_of_tails
-          DO  n = number_of_particles+1, number_of_particles+trnp_count_recv
-             IF ( particles(n)%tail_id /= 0 )  THEN
-                nn = nn + 1
-                particles(n)%tail_id = nn
-             ENDIF
-          ENDDO
-
-       ENDIF
-
-!       number_of_particles = number_of_particles + trnp_count_recv
-!       number_of_tails     = number_of_tails     + trnpt_count_recv
-
 !
 !--    Send back boundary, receive front boundary
@@ -876 +644 @@
        DEALLOCATE(rvsp)
 
-       IF ( use_particle_tails )  THEN
-
-          CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pnorth, 0, &
-                             trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
-                             comm2d, status, ierr )
-
-          IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
-          THEN
-             IF ( netcdf_data_format < 3 )  THEN
-                message_string = 'maximum_number_of_tails ' //   &
-                                 'needs to be increased ' //     &
-                                 '&but this is not allowed wi'// &
-                                 'th NetCDF output switched on'
-                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
-             ELSE
-                CALL lpm_extend_tail_array( trspt_count_recv )
-             ENDIF
-          ENDIF
-
-          CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL,      &
-                             pnorth, 1,                                        &
-                             particle_tail_coordinates(1,1,number_of_tails+1), &
-                             trspt_count_recv*tlength, MPI_REAL, psouth, 1,    &
-                             comm2d, status, ierr )
-!
-!--       Update the tail ids for the transferred particles
-          nn = number_of_tails
-          DO  n = number_of_particles+1, number_of_particles+trsp_count_recv
-             IF ( particles(n)%tail_id /= 0 )  THEN
-                nn = nn + 1
-                particles(n)%tail_id = nn
-             ENDIF
-          ENDDO
-
-       ENDIF
-
        number_of_particles = number_of_particles + trsp_count_recv
-       number_of_tails     = number_of_tails     + trspt_count_recv
-
-       IF ( use_particle_tails )  THEN
-          DEALLOCATE( trspt, trnpt )
-       ENDIF
+
        DEALLOCATE( trsp, trnp )
 
@@ -928 +656 @@
     DO  n = 1, number_of_particles
 
-       nn = particles(n)%tail_id
-
        IF ( particles(n)%x < -0.5_wp * dx )  THEN
 
@@ -936 +662 @@
 !--          Cyclic boundary. Relevant coordinate has to be changed.
              particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                i = particles(n)%tailpoints
-                particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
-                                             particle_tail_coordinates(1:i,1,nn)
-             ENDIF
+
           ELSEIF ( ibc_par_lr == 1 )  THEN
 !
@@ -946 +668 @@
              particles(n)%particle_mask = .FALSE.
              deleted_particles = deleted_particles + 1
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                tail_mask(nn) = .FALSE.
-                deleted_tails = deleted_tails + 1
-             ENDIF
+
           ELSEIF ( ibc_par_lr == 2 )  THEN
 !
@@ -963 +682 @@
 !--          Cyclic boundary. Relevant coordinate has to be changed.
              particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                i = particles(n)%tailpoints
-                particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
-                                             particle_tail_coordinates(1:i,1,nn)
-             ENDIF
+
           ELSEIF ( ibc_par_lr == 1 )  THEN
 !
@@ -973 +688 @@
              particles(n)%particle_mask = .FALSE.
              deleted_particles = deleted_particles + 1
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                tail_mask(nn) = .FALSE.
-                deleted_tails = deleted_tails + 1
-             ENDIF
+
           ELSEIF ( ibc_par_lr == 2 )  THEN
 !
@@ -992 +704 @@
 !--          Cyclic boundary. Relevant coordinate has to be changed.
              particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                i = particles(n)%tailpoints
-                particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
-                                             particle_tail_coordinates(1:i,2,nn)
-             ENDIF
+
           ELSEIF ( ibc_par_ns == 1 )  THEN
 !
@@ -1002 +710 @@
              particles(n)%particle_mask = .FALSE.
              deleted_particles = deleted_particles + 1
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                tail_mask(nn) = .FALSE.
-                deleted_tails = deleted_tails + 1
-             ENDIF
+
           ELSEIF ( ibc_par_ns == 2 )  THEN
 !
@@ -1019 +724 @@
 !--          Cyclic boundary. Relevant coordinate has to be changed.
              particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                i = particles(n)%tailpoints
-                particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
-                                             particle_tail_coordinates(1:i,2,nn)
-             ENDIF
+
           ELSEIF ( ibc_par_ns == 1 )  THEN
 !
@@ -1029 +730 @@
              particles(n)%particle_mask = .FALSE.
              deleted_particles = deleted_particles + 1
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                tail_mask(nn) = .FALSE.
-                deleted_tails = deleted_tails + 1
-             ENDIF
+
           ELSEIF ( ibc_par_ns == 2 )  THEN
 !
@@ -1287 +985 @@
     INTEGER(iwp), INTENT(in)                       ::  j              !<
     INTEGER(iwp), INTENT(in)                       ::  k              !<
-    INTEGER(iwp), INTENT(in), optional             ::  size_in        !<
+    INTEGER(iwp), INTENT(in), OPTIONAL             ::  size_in        !<
 
     INTEGER(iwp)                                   :: old_size        !<
@@ -1300 +998 @@
        new_size = size_in
     ELSE
-       new_size = old_size * ( 1.0 + alloc_factor / 100.0 )
+       new_size = old_size * ( 1.0_wp + alloc_factor / 100.0_wp )
     ENDIF
 

palm/trunk/SOURCE/lpm_init.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Unused variables removed.
 !
 ! Former revisions:
@@ -113 +113 @@
                dz, initializing_actions, message_string, ocean, simulated_time
 
-    USE dvrp_variables,                                                        &
-        ONLY:  particle_color, particle_dvrpsize
-
     USE grid_variables,                                                        &
         ONLY:  ddx, dx, ddy, dy
@@ -134 +131 @@
         ONLY:   alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,        &
                 block_offset, block_offset_def, collision_kernel,              &
-                density_ratio, dvrp_psize, grid_particles,                     &
+                density_ratio, grid_particles,                                 &
                 initial_weighting_factor, ibc_par_b, ibc_par_lr, ibc_par_ns,   &
                 ibc_par_t, iran_part, log_z_z0,                                &
                 max_number_of_particle_groups, maximum_number_of_particles,    &
-                maximum_number_of_tailpoints, maximum_number_of_tails,         &
-                minimum_tailpoint_distance, min_nr_particle,                   &
-                mpi_particle_type, new_tail_id,                                &
-                number_of_initial_tails, number_of_particles,                  &
+                min_nr_particle, mpi_particle_type,                            &
+                number_of_particles,                                           &
                 number_of_particle_groups, number_of_sublayers,                &
-                number_of_tails, offset_ocean_nzt, offset_ocean_nzt_m1,        &
+                offset_ocean_nzt, offset_ocean_nzt_m1,                         &
                 particles, particle_advection_start, particle_groups,          &
                 particle_groups_type, particles_per_point,                     &
-                particle_tail_coordinates,  particle_type, pdx, pdy, pdz,      &
+                particle_type, pdx, pdy, pdz,                                  &
                 prt_count, psb, psl, psn, psr, pss, pst,                       &
                 radius, random_start_position, read_particles_from_restartfile,&
-                seed_follows_topography, skip_particles_for_tail, sort_count,  &
-                tail_mask, total_number_of_particles, total_number_of_tails,   &
-                use_particle_tails, use_sgs_for_particles,                     &
+                seed_follows_topography, sort_count,                           &
+                total_number_of_particles,                                     &
+                use_sgs_for_particles,                                         &
                 write_particle_statistics, uniform_particles, zero_particle,   &
                 z0_av_global
@@ -192 +187 @@
 
     INTEGER(iwp) ::  i                           !<
-    INTEGER(iwp) ::  ip                          !<
     INTEGER(iwp) ::  j                           !<
-    INTEGER(iwp) ::  jp                          !<
     INTEGER(iwp) ::  k                           !<
-    INTEGER(iwp) ::  kp                          !<
-    INTEGER(iwp) ::  n                           !<
-    INTEGER(iwp) ::  nn                          !<
 
 #if defined( __parallel )
@@ -205 +195 @@
     INTEGER(iwp), DIMENSION(3) ::  types         !<
 #endif
-    LOGICAL ::  uniform_particles_l              !<
 
     REAL(wp) ::  height_int                      !<
@@ -216 +205 @@
 !-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
 !-- particle_attributes). Integer length is 4 byte, Real is 8 byte
-#if defined( __twocachelines )
-    blocklengths(1)  =  7;  blocklengths(2)  =  18;  blocklengths(3)  =   1
-    displacements(1) =  0;  displacements(2) =  64;  displacements(3) = 128
-
-    types(1) = MPI_REAL                               ! 64 bit words
-    types(2) = MPI_INTEGER                            ! 32 Bit words
-    types(3) = MPI_UB
-#else
     blocklengths(1)  = 19;  blocklengths(2)  =   6;  blocklengths(3)  =   1
     displacements(1) =  0;  displacements(2) = 152;  displacements(3) = 176
@@ -230 +211 @@
     types(2) = MPI_INTEGER
     types(3) = MPI_UB
-#endif
     CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, &
                           mpi_particle_type, ierr )
@@ -293 +273 @@
 !
 !-- Allocate arrays required for calculating particle SGS velocities
-    IF ( use_sgs_for_particles )  THEN
+    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
        ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
                  de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
@@ -429 +409 @@
 !--    Allocate particle arrays and set attributes of the initial set of
 !--    particles, which can be also periodically released at later times.
-!--    Also allocate array for particle tail coordinates, if needed.
        ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
                  grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )
@@ -443 +422 @@
 !--    occur within restart runs). The reason for this is still not clear
 !--    and may be presumably caused by errors in the respective user-interface.
-#if defined( __twocachelines )
-       zero_particle = particle_type( 0.0_wp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp,  &
-                                      0.0_sp, 0.0_sp, 0.0_wp, 0.0_wp, 0.0_wp,  &
-                                      0, .FALSE., 0.0_wp, 0.0_wp, 0.0_wp,      &
-                                      0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp,  &
-                                      0.0_sp, 0, 0, 0, -1)
-#else
        zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
@@ -455 +427 @@
                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, &
                                       0, .FALSE., -1)
-#endif
+
        particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
-
-!
-!--    Set the default particle size used for dvrp plots
-       IF ( dvrp_psize == 9999999.9_wp )  dvrp_psize = 0.2_wp * dx
 
 !
@@ -506 +474 @@
        ENDIF
 
-!
-!--    Check if particles are really uniform in color and radius (dvrp_size)
-!--    (uniform_particles is preset TRUE)
-       IF ( uniform_particles )  THEN
-          DO  ip = nxl, nxr
-             DO  jp = nys, nyn
-                DO  kp = nzb+1, nzt
-
-                   n = prt_count(kp,jp,ip)
-                   IF ( MINVAL( grid_particles(kp,jp,ip)%particles(1:n)%dvrp_psize  ) ==     &
-                        MAXVAL( grid_particles(kp,jp,ip)%particles(1:n)%dvrp_psize  )  .AND. &
-                        MINVAL( grid_particles(kp,jp,ip)%particles(1:n)%class ) ==           &
-                        MAXVAL( grid_particles(kp,jp,ip)%particles(1:n)%class ) )  THEN
-                      uniform_particles_l = .TRUE.
-                   ELSE
-                      uniform_particles_l = .FALSE.
-                   ENDIF
-
-                ENDDO
-             ENDDO
-          ENDDO
-
-#if defined( __parallel )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( uniform_particles_l, uniform_particles, 1,       &
-                              MPI_LOGICAL, MPI_LAND, comm2d, ierr )
-#else
-          uniform_particles = uniform_particles_l
-#endif
-
-       ENDIF
-
-!
-!--    Particles will probably become none-uniform, if their size and color
-!--    will be determined by flow variables
-       IF ( particle_color /= 'none'  .OR.  particle_dvrpsize /= 'none' )  THEN
-          uniform_particles = .FALSE.
-       ENDIF
-
-! !kk    Not implemented aft individual particle array fort every gridcell
-! !
-! !--    Set the beginning of the particle tails and their age
-!        IF ( use_particle_tails )  THEN
-! !
-! !--       Choose the maximum number of tails with respect to the maximum number
-! !--       of particles and skip_particles_for_tail
-!           maximum_number_of_tails = maximum_number_of_particles / &
-!                                     skip_particles_for_tail
-! 
-! !
-! !--       Create a minimum number of tails in case that there is no tail
-! !--       initially (otherwise, index errors will occur when adressing the
-! !--       arrays below)
-!           IF ( maximum_number_of_tails == 0 )  maximum_number_of_tails = 10
-! 
-!           ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &
-!                     maximum_number_of_tails),                                 &
-!                     new_tail_id(maximum_number_of_tails),                     &
-!                     tail_mask(maximum_number_of_tails) )
-! 
-!           particle_tail_coordinates  = 0.0_wp
-!           minimum_tailpoint_distance = minimum_tailpoint_distance**2
-!           number_of_initial_tails    = number_of_tails
-! 
-!           nn = 0
-!           DO  n = 1, number_of_particles
-! !
-! !--          Only for those particles marked above with a provisional tail_id
-! !--          tails will be created. Particles now get their final tail_id.
-!              IF ( particles(n)%tail_id /= 0 )  THEN
-! 
-!                 nn = nn + 1
-!                 particles(n)%tail_id = nn
-! 
-!                 particle_tail_coordinates(1,1,nn) = particles(n)%x
-!                 particle_tail_coordinates(1,2,nn) = particles(n)%y
-!                 particle_tail_coordinates(1,3,nn) = particles(n)%z
-!                 particle_tail_coordinates(1,4,nn) = particles(n)%class
-!                 particles(n)%tailpoints = 1
-!                 IF ( minimum_tailpoint_distance /= 0.0_wp )  THEN
-!                    particle_tail_coordinates(2,1,nn) = particles(n)%x
-!                    particle_tail_coordinates(2,2,nn) = particles(n)%y
-!                    particle_tail_coordinates(2,3,nn) = particles(n)%z
-!                    particle_tail_coordinates(2,4,nn) = particles(n)%class
-!                    particle_tail_coordinates(1:2,5,nn) = 0.0_wp
-!                    particles(n)%tailpoints = 2
-!                 ENDIF
-! 
-!              ENDIF
-!           ENDDO
-!        ENDIF
-! 
-! !
-! !--    Plot initial positions of particles (only if particle advection is
-! !--    switched on from the beginning of the simulation (t=0))
-!        IF ( particle_advection_start == 0.0_wp )  CALL data_output_dvrp
-
     ENDIF
 
@@ -641 +512 @@
     INTEGER(iwp)               ::  n           !<
     INTEGER(iwp)               ::  new_size    !<
-    INTEGER(iwp)               ::  nn          !<
 
     INTEGER(iwp), INTENT(IN)   ::  phase       !<
@@ -691 +561 @@
 
                             n = n + 1
-#if defined( __twocachelines )
-                            tmp_particle%x             = pos_x
-                            tmp_particle%y             = pos_y
-                            tmp_particle%z             = pos_z
-                            tmp_particle%age           = 0.0_sp
-                            tmp_particle%age_m         = 0.0_sp
-                            tmp_particle%dt_sum        = 0.0_wp
-                            tmp_particle%dvrp_psize    = dvrp_psize
-                            tmp_particle%e_m           = 0.0_sp
-                            IF ( curvature_solution_effects )  THEN
-!
-!--                            Initial values (internal timesteps, derivative)
-!--                            for Rosenbrock method
-                               tmp_particle%rvar1      = 1.0E-12_wp
-                               tmp_particle%rvar2      = 1.0E-3_wp
-                               tmp_particle%rvar3      = -9999999.9_wp
-                            ELSE
-!
-!--                            Initial values for SGS velocities
-                               tmp_particle%rvar1      = 0.0_wp
-                               tmp_particle%rvar2      = 0.0_wp
-                               tmp_particle%rvar3      = 0.0_wp
-                            ENDIF
-                            tmp_particle%speed_x       = 0.0_sp
-                            tmp_particle%speed_y       = 0.0_sp
-                            tmp_particle%speed_z       = 0.0_sp
-                            tmp_particle%origin_x      = pos_x
-                            tmp_particle%origin_y      = pos_y
-                            tmp_particle%origin_z      = pos_z
-                            tmp_particle%radius        = particle_groups(i)%radius
-                            tmp_particle%weight_factor = initial_weighting_factor
-                            tmp_particle%class         = 1
-                            tmp_particle%group         = i
-                            tmp_particle%tailpoints    = 0
-                            tmp_particle%particle_mask = .TRUE.
-#else
                             tmp_particle%x             = pos_x
                             tmp_particle%y             = pos_y
@@ -733 +567 @@
                             tmp_particle%age_m         = 0.0_wp
                             tmp_particle%dt_sum        = 0.0_wp
-                            tmp_particle%dvrp_psize    = dvrp_psize
+                            tmp_particle%dvrp_psize    = 0.0_wp !unused
                             tmp_particle%e_m           = 0.0_wp
                             IF ( curvature_solution_effects )  THEN
@@ -759 +593 @@
                             tmp_particle%class         = 1
                             tmp_particle%group         = i
-                            tmp_particle%tailpoints    = 0
+                            tmp_particle%tailpoints    = 0     !unused
                             tmp_particle%particle_mask = .TRUE.
-#endif
-                            IF ( use_particle_tails  .AND. &
-                                 MOD( n, skip_particles_for_tail ) == 0 )  THEN
-                               number_of_tails         = number_of_tails + 1
-!
-!--                            This is a temporary provisional setting (see
-!--                            further below!)
-                               tmp_particle%tail_id    = number_of_tails
-                            ELSE
-                               tmp_particle%tail_id    = 0
-                            ENDIF
+                            tmp_particle%tail_id       = 0     !unused
+
 !
 !--                         Determine the grid indices of the particle position
@@ -913 +738 @@
     ENDDO
 !
-!-- Calculate the number of particles and tails of the total domain
+!-- Calculate the number of particles of the total domain
 #if defined( __parallel )
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &
     MPI_INTEGER, MPI_SUM, comm2d, ierr )
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( number_of_tails, total_number_of_tails, 1, &
-    MPI_INTEGER, MPI_SUM, comm2d, ierr )
 #else
     total_number_of_particles = number_of_particles
-    total_number_of_tails     = number_of_tails
 #endif
 

palm/trunk/SOURCE/lpm_init_sgs_tke.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Unused variables removed.
 ! 
 ! Former revisions:
@@ -63 +63 @@
 
     USE indices,                                                               &
-        ONLY:  nbgp, ngp_2dh_outer, nx, nxl, nxr, ny, nyn, nys, nz, nzb,       &
-                      nzb_s_inner, nzb_s_outer, nzt
+        ONLY:  nbgp, ngp_2dh_outer, nxl, nxr, nyn, nys, nzb, nzb_s_inner,      &
+               nzb_s_outer, nzt
 
     USE kinds

palm/trunk/SOURCE/lpm_pack_arrays.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Tails removed. Unused variables removed.
 !
 ! Former revisions:
@@ -51 +51 @@
 ! Description:
 ! ------------
-!> Pack particle and tail arrays, which means eliminate those elements marked for
+!> Pack particle arrays, which means eliminate those elements marked for
 !> deletion and move data with higher index values to these free indices.
 !> Determine the new number of particles.
@@ -59 +59 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  deleted_tails, grid_particles, new_tail_id,                     &
-               number_of_particles, number_of_tails, offset_ocean_nzt,         &
-               particles, particle_type, prt_count,                            &
-               particle_tail_coordinates, tail_mask, use_particle_tails
+        ONLY:  grid_particles, number_of_particles, offset_ocean_nzt,          &
+               particles, particle_type, prt_count
 
     PRIVATE
@@ -127 +125 @@
 
        INTEGER(iwp) ::  n       !<
-       INTEGER(iwp) ::  nd      !<
        INTEGER(iwp) ::  nn      !<
 !
@@ -158 +155 @@
        number_of_particles = nn
 
-!
-!-- particle tails are currently not available
-!
-!-- Handle tail array in the same way, store the new tail ids and re-assign it
-!-- to the respective particles
-!    IF ( use_particle_tails )  THEN
-!
-!       nn = 0
-!       nd = 0
-!
-!       DO  n = 1, number_of_tails
-!
-!          IF ( tail_mask(n) )  THEN
-!             nn = nn + 1
-!             particle_tail_coordinates(:,:,nn) = &
-!                                                particle_tail_coordinates(:,:,n)
-!             new_tail_id(n) = nn
-!          ELSE
-!             nd = nd + 1
-!          ENDIF
-!
-!       ENDDO
-!
-!       DO  n = 1, number_of_particles
-!          IF ( particles(n)%tail_id /= 0 )  THEN
-!             particles(n)%tail_id = new_tail_id(particles(n)%tail_id)
-!          ENDIF
-!       ENDDO
-!
-!    ENDIF
-
-!
-!-- The number of deleted tails has been determined in routines
-!-- lpm_boundary_conds and lpm_exchange_horiz
-!   number_of_tails = number_of_tails - deleted_tails
-
-
     END SUBROUTINE lpm_pack_arrays
 
@@ -206 +166 @@
       USE control_parameters,                                                  &
           ONLY: dz
-
-      USE indices,                                                             &
-          ONLY: nxl, nxr, nys, nyn, nzb, nzt
 
       USE kinds

palm/trunk/SOURCE/lpm_read_restart_file.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Tails removed. Unused variables removed.
 ! 
 ! Former revisions:
@@ -54 +54 @@
 
     USE indices,                                                               &
-        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, nzb, nzt
+        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
 
     USE kinds
@@ -64 +64 @@
         ONLY:  alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,         &
                grid_particles, maximum_number_of_particles,                    &
-               maximum_number_of_tailpoints, maximum_number_of_tails,          &
-               min_nr_particle, new_tail_id, number_of_particles,              &
-               number_of_particle_groups, number_of_tails, particles,          &
-               particle_groups, particle_tail_coordinates, particle_type,      &
-               prt_count, sort_count, tail_mask, time_prel,                    &
-               time_write_particle_data, uniform_particles,                    &
-               use_particle_tails, zero_particle
+               min_nr_particle, number_of_particles, number_of_particle_groups,&
+               particle_groups, particle_type, prt_count, time_prel,           &
+               time_write_particle_data, uniform_particles, zero_particle
 
     USE pegrid
@@ -100 +96 @@
 !-- First compare the version numbers
     READ ( 90 )  version_on_file
-    particle_binary_version = '3.2'
+    particle_binary_version = '4.0'
     IF ( TRIM( version_on_file ) /= TRIM( particle_binary_version ) )  THEN
        message_string = 'version mismatch concerning data from prior ' // &
@@ -114 +110 @@
 !-- min_nr_particle, the remainder is initialized by zero_particle to avoid
 !-- errors.
-#if defined( __twocachelines )
-    zero_particle = particle_type( 0.0_wp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp,  &
-                                   0.0_sp, 0.0_sp, 0.0_wp, 0.0_wp, 0.0_wp,  &
-                                   0, .FALSE., 0.0_wp, 0.0_wp, 0.0_wp,      &
-                                   0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp,  &
-                                   0.0_sp, 0, 0, 0, -1)
-#else
     zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
                                    0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
@@ -126 +115 @@
                                    0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, &
                                    0, .FALSE., -1)
-#endif
-
 !
 !-- Read some particle parameters and the size of the particle arrays,
 !-- allocate them and read their contents.
     READ ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,                  &
-                 maximum_number_of_tailpoints, maximum_number_of_tails,     &
-                 number_of_particle_groups, number_of_tails,                &
-                 particle_groups, time_prel, time_write_particle_data,      &
-                 uniform_particles
+                 number_of_particle_groups, particle_groups, time_prel,     &
+                 time_write_particle_data, uniform_particles
 
     ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                     &
@@ -177 +162 @@
     ENDDO
 
-!
-!-- particle tails currently not available
-!    IF ( use_particle_tails )  THEN
-!       ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &
-!                 maximum_number_of_tails),                                 &
-!                 new_tail_id(maximum_number_of_tails),                     &
-!                 tail_mask(maximum_number_of_tails) )
-!       READ ( 90 )  particle_tail_coordinates
-!    ENDIF
-
     CLOSE ( 90 )
 !

palm/trunk/SOURCE/lpm_set_attributes.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Unused variables removed.
 ! 
 ! Former revisions:
@@ -67 +67 @@
 
     USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, u_gtrans, v_gtrans, dz
+        ONLY:  u_gtrans, v_gtrans, dz
 
     USE cpulog,                                                                &
@@ -77 +77 @@
 
     USE grid_variables,                                                        &
-        ONLY:  ddx, dx, ddy, dy
+        ONLY:  dx, dy
 
     USE indices,                                                               &
@@ -85 +85 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  block_offset, grid_particles, number_of_particles,              &
-               offset_ocean_nzt, particles, prt_count
+        ONLY:  block_offset, grid_particles, number_of_particles, particles,   &
+               prt_count
 
     USE pegrid

palm/trunk/SOURCE/lpm_write_exchange_statistics.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Unused variables removed.
 ! 
 ! Former revisions:
@@ -62 +62 @@
     USE particle_attributes,                                                   &
         ONLY:  grid_particles, maximum_number_of_particles,                    &
-               number_of_particles, particles, prt_count,                      &
+               number_of_particles, prt_count,                                 &
                trlp_count_sum, trlp_count_recv_sum, trnp_count_sum,            &
                trnp_count_recv_sum, trrp_count_sum, trrp_count_recv_sum,       &

palm/trunk/SOURCE/lpm_write_restart_file.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Tails removed. Unused variables removed.
 ! 
 ! Former revisions:
@@ -50 +50 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lpm_write_restart_file
- 
 
-    USE control_parameters,                                                    &
-        ONLY:  io_blocks, io_group
 
     USE indices,                                                               &
@@ -62 +59 @@
     USE particle_attributes,                                                   &
         ONLY:  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, grid_particles,       &
-               maximum_number_of_tails, maximum_number_of_tailpoints,          &
                number_of_particles, number_of_particle_groups,                 &
-               number_of_tails, particles, particle_groups,                    &
-               particle_tail_coordinates, prt_count, time_prel,                &
-               time_write_particle_data, uniform_particles,                    &
-               use_particle_tails, zero_particle
+               particles, particle_groups, prt_count, time_prel,               &
+               time_write_particle_data, uniform_particles
 
     USE pegrid
@@ -75 +69 @@
     CHARACTER (LEN=10) ::  particle_binary_version   !<
 
-    INTEGER(iwp) ::  i                               !<
     INTEGER(iwp) ::  ip                              !<
     INTEGER(iwp) ::  jp                              !<
@@ -97 +90 @@
     ENDIF
 
-!    DO  i = 0, io_blocks-1
-
-!       IF ( i == io_group )  THEN
+!
+!-- Write the version number of the binary format.
+!-- Attention: After changes to the following output commands the version
+!-- ---------  number of the variable particle_binary_version must be
+!--            changed! Also, the version number and the list of arrays
+!--            to be read in lpm_read_restart_file must be adjusted
+!--            accordingly.
+    particle_binary_version = '4.0'
+    WRITE ( 90 )  particle_binary_version
 
 !
-!--       Write the version number of the binary format.
-!--       Attention: After changes to the following output commands the version
-!--       ---------  number of the variable particle_binary_version must be
-!--                  changed! Also, the version number and the list of arrays
-!--                  to be read in lpm_read_restart_file must be adjusted
-!--                  accordingly.
-          particle_binary_version = '3.2'
-          WRITE ( 90 )  particle_binary_version
+!-- Write some particle parameters, the size of the particle arrays as
+!-- well as other dvrp-plot variables.
+    WRITE ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,              &
+                  number_of_particle_groups, particle_groups, time_prel, &
+                  time_write_particle_data, uniform_particles
 
-!
-!--       Write some particle parameters, the size of the particle arrays as
-!--       well as other dvrp-plot variables.
-          WRITE ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,              &
-                        maximum_number_of_tailpoints, maximum_number_of_tails, &
-                        number_of_particle_groups, number_of_tails,            &
-                        particle_groups, time_prel, time_write_particle_data,  &
-                        uniform_particles
+    WRITE ( 90 )  prt_count
+          
+    DO  ip = nxl, nxr
+       DO  jp = nys, nyn
+          DO  kp = nzb+1, nzt
+             number_of_particles = prt_count(kp,jp,ip)
+             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+             IF ( number_of_particles <= 0 )  CYCLE
+             WRITE ( 90 )  particles
+          ENDDO
+       ENDDO
+    ENDDO
 
-          WRITE ( 90 )  prt_count
-          
-          DO  ip = nxl, nxr
-             DO  jp = nys, nyn
-                DO  kp = nzb+1, nzt
-                   number_of_particles = prt_count(kp,jp,ip)
-                   particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
-                   IF ( number_of_particles <= 0 )  CYCLE
-                   WRITE ( 90 )  particles
-                ENDDO
-             ENDDO
-          ENDDO
-
-!
-!--       particle tails currently not available
-!          IF ( use_particle_tails )  THEN
-!             WRITE ( 90 )  particle_tail_coordinates
-!          ENDIF
-
-          CLOSE ( 90 )
-
-!       ENDIF
+    CLOSE ( 90 )
 
 #if defined( __parallel )
@@ -147 +126 @@
 #endif
 
- !   ENDDO
-
 
  END SUBROUTINE lpm_write_restart_file

palm/trunk/SOURCE/microphysics.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Unused variables removed. 
+!
+! Kessler scheme integrated.
 ! 
 ! Former revisions:
@@ -109 +111 @@
     END INTERFACE autoconversion
 
+    INTERFACE autoconversion_kessler
+       MODULE PROCEDURE autoconversion_kessler
+       MODULE PROCEDURE autoconversion_kessler_ij
+    END INTERFACE autoconversion_kessler
+
     INTERFACE accretion
        MODULE PROCEDURE accretion
@@ -151 +158 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  hyp, nr, pt, pt_init, q, qc, qr, zu
+           ONLY:  hyp, pt_init, zu
 
        USE cloud_parameters,                                                   &
-           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
+           ONLY:  cp, hyrho, prr, pt_d_t, r_d, t_d_pt
 
        USE control_parameters,                                                 &
            ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
-                  g, intermediate_timestep_count,                              &
-                  large_scale_forcing, lsf_surf, precipitation, pt_surface,    &
-                  rho_surface,surface_pressure
+                  g, intermediate_timestep_count, large_scale_forcing,         &
+                  lsf_surf, microphysics_kessler, microphysics_seifert,        &
+                  pt_surface, rho_surface,surface_pressure
 
        USE indices,                                                            &
@@ -172 +179 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
        INTEGER(iwp) ::  k                 !<
 
@@ -193 +198 @@
              hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
           ENDDO
+
 !
 !--       Compute reference density
@@ -207 +213 @@
 
 !
+!--    Reset precipitation rate
+       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
+
+!
 !--    Compute cloud physics
-       IF ( precipitation )  THEN
+       IF ( microphysics_kessler )  THEN
+
+          CALL autoconversion_kessler
+
+       ELSEIF ( microphysics_seifert )  THEN
+
           CALL adjust_cloud
           CALL autoconversion
@@ -215 +230 @@
           CALL evaporation_rain
           CALL sedimentation_rain
+
+          IF ( drizzle )  CALL sedimentation_cloud
+
        ENDIF
 
-       IF ( drizzle )  CALL sedimentation_cloud
-
-       IF ( precipitation )  THEN
-          CALL calc_precipitation_amount
-       ENDIF
+       CALL calc_precipitation_amount
 
     END SUBROUTINE microphysics_control
@@ -238 +252 @@
 
        USE cloud_parameters,                                                   &
-           ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
+           ONLY:  eps_sb, xrmin, xrmax, hyrho
 
        USE cpulog,                                                             &
@@ -244 +258 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
 
        USE kinds
@@ -303 +317 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
 
        USE kinds
@@ -323 +337 @@
        REAL(wp)     ::  rc                !<
        REAL(wp)     ::  re_lambda         !<
-       REAL(wp)     ::  selfcoll          !<
        REAL(wp)     ::  sigma_cc          !<
        REAL(wp)     ::  tau_cloud         !<
@@ -417 +430 @@
 ! Description:
 ! ------------
+!> Autoconversion process (Kessler, 1969).
+!------------------------------------------------------------------------------!
+    SUBROUTINE autoconversion_kessler
+
+       USE arrays_3d,                                                          &
+           ONLY:  dzw, pt, q, qc
+
+       USE cloud_parameters,                                                   &
+           ONLY:  l_d_cp, pt_d_t, prec_time_const, prr, ql_crit
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_micro
+
+       USE indices,                                                            &
+           ONLY:  nxl, nxr, nyn, nys, nzb_2d, nzt
+
+       USE kinds
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !<
+       INTEGER(iwp) ::  j !<
+       INTEGER(iwp) ::  k !<
+
+       REAL(wp)    ::  dqdt_precip !<
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_2d(j,i)+1, nzt
+
+                IF ( qc(k,j,i) > ql_crit )  THEN
+                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
+                ELSE
+                   dqdt_precip = 0.0_wp
+                ENDIF
+
+                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro
+                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro
+                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k)
+
+!
+!--             Compute the rain rate
+                prr(nzb_2d(j,i)+1,j,i) = prr(nzb_2d(j,i)+1,j,i) + dqdt_precip * dzw(k)
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+   END SUBROUTINE autoconversion_kessler
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Accretion rate (Seifert and Beheng, 2006).
 !------------------------------------------------------------------------------!
@@ -434 +502 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
 
        USE kinds
@@ -448 +516 @@
        REAL(wp)     ::  phi_ac            !<
        REAL(wp)     ::  tau_cloud         !<
-       REAL(wp)     ::  xc                !<
 
        CALL cpu_log( log_point_s(56), 'accretion', 'start' )
@@ -516 +583 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
 
        USE kinds
@@ -579 +646 @@
 
        USE cloud_parameters,                                                   &
-           ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
-                  dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,   &
-                  l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,      &
-                  t_d_pt, ventilation_effect
+           ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term,              &
+                  diff_coeff_l, dpirho_l, eps_sb, hyrho, kin_vis_air,          &
+                  l_d_cp, l_d_r, l_v, r_v, schmidt_p_1d3,                      &
+                  thermal_conductivity_l, t_d_pt, ventilation_effect
 
        USE constants,                                                          &
@@ -594 +661 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
 
        USE kinds
@@ -742 +809 @@
            ONLY:  eps_sb, hyrho, l_d_cp, nc_const, prr, pt_d_t, sed_qc_const
 
-       USE constants,                                                          &
-           ONLY:  pi
-
        USE control_parameters,                                                 &
            ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
-                  intermediate_timestep_count, precipitation
+                  intermediate_timestep_count
 
        USE cpulog,                                                             &
@@ -798 +862 @@
 !
 !--             Compute the precipitation rate due to cloud (fog) droplets
-                IF ( precipitation )  THEN
-                   IF ( call_microphysics_at_all_substeps )  THEN
-                      prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)          &
-                                   * weight_substep(intermediate_timestep_count)
-                   ELSE
-                      prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
-                   ENDIF
+                IF ( call_microphysics_at_all_substeps )  THEN
+                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
+                                * weight_substep(intermediate_timestep_count)
+                ELSE
+                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
                 ENDIF
 
@@ -828 +890 @@
 
        USE cloud_parameters,                                                   &
-           ONLY:  a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,        &
-                  limiter_sedimentation, l_d_cp, prr,  pt_d_t, stp
+           ONLY:  a_term, b_term, c_term, dpirho_l, eps_sb, hyrho,             &
+                  limiter_sedimentation, l_d_cp, prr, pt_d_t
 
        USE control_parameters,                                                 &
@@ -857 +919 @@
        REAL(wp)     ::  d_min                      !<
        REAL(wp)     ::  dr                         !<
-       REAL(wp)     ::  dt_sedi                    !<
        REAL(wp)     ::  flux                       !<
        REAL(wp)     ::  lambda_r                   !<
@@ -865 +926 @@
        REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
        REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr     !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_qr     !<
        REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
        REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
@@ -876 +935 @@
        CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
 
-       IF ( intermediate_timestep_count == 1 )  prr(:,:,:) = 0.0_wp
 !
 !--    Compute velocities 
@@ -1105 +1163 @@
 
        USE cloud_parameters,                                                   &
-           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
+           ONLY:  cp, hyrho, nc_const, prr, pt_d_t, r_d, t_d_pt
 
        USE control_parameters,                                                 &
            ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
                   g, intermediate_timestep_count, large_scale_forcing,         &
-                  lsf_surf, precipitation, pt_surface,                         &
-                  rho_surface,surface_pressure
+                  lsf_surf, microphysics_seifert, microphysics_kessler,        &
+                  pt_surface, rho_surface, surface_pressure
 
        USE indices,                                                            &
@@ -1162 +1220 @@
        qc_1d(:) = qc(:,j,i)
        nc_1d(:) = nc_const
-       IF ( precipitation )  THEN
+       IF ( microphysics_seifert )  THEN
           qr_1d(:) = qr(:,j,i)
           nr_1d(:) = nr(:,j,i)
@@ -1168 +1226 @@
 
 !
+!--    Reset precipitation rate
+       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
+
+!
 !--    Compute cloud physics
-       IF ( precipitation )  THEN
-          CALL adjust_cloud( i,j ) 
+       IF( microphysics_kessler )  THEN
+
+          CALL autoconversion_kessler( i,j )
+
+       ELSEIF ( microphysics_seifert )  THEN
+
+          CALL adjust_cloud( i,j )
           CALL autoconversion( i,j )
           CALL accretion( i,j )
@@ -1176 +1243 @@
           CALL evaporation_rain( i,j )
           CALL sedimentation_rain( i,j )
+
+          IF ( drizzle )  CALL sedimentation_cloud( i,j )
+
        ENDIF
 
-       IF ( drizzle )  CALL sedimentation_cloud( i,j )
-
-       IF ( precipitation )  THEN
-          CALL calc_precipitation_amount( i,j )
-       ENDIF
+       CALL calc_precipitation_amount( i,j )
 
 !
@@ -1188 +1254 @@
        q(:,j,i)  = q_1d(:)
        pt(:,j,i) = pt_1d(:)
-       IF ( precipitation )  THEN
+       IF ( microphysics_seifert )  THEN
           qr(:,j,i) = qr_1d(:)
           nr(:,j,i) = nr_1d(:)
@@ -1210 +1276 @@
 
        USE cloud_parameters,                                                   &
-           ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
+           ONLY:  eps_sb, xrmin, xrmax, hyrho
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nzb_s_inner, nzt
 
        USE kinds
@@ -1258 +1324 @@
        USE cloud_parameters,                                                   &
            ONLY:  a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3,        &
-                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
+                  c_const, dpirho_l, eps_sb, hyrho, kin_vis_air, k_cc, x0
 
        USE control_parameters,                                                 &
@@ -1267 +1333 @@
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nzb_s_inner, nzt
 
        USE kinds
@@ -1287 +1353 @@
        REAL(wp)     ::  rc                !<
        REAL(wp)     ::  re_lambda         !<
-       REAL(wp)     ::  selfcoll          !<
        REAL(wp)     ::  sigma_cc          !<
        REAL(wp)     ::  tau_cloud         !<
@@ -1366 +1431 @@
     END SUBROUTINE autoconversion_ij
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Autoconversion process (Kessler, 1969).
+!------------------------------------------------------------------------------!
+    SUBROUTINE autoconversion_kessler_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  dzw, pt_1d, q_1d, qc_1d
+
+       USE cloud_parameters,                                                   &
+           ONLY:  l_d_cp, pt_d_t, prec_time_const, prr, ql_crit
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_micro
+
+       USE indices,                                                            &
+           ONLY:  nzb_2d, nzt
+
+       USE kinds
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !<
+       INTEGER(iwp) ::  j !<
+       INTEGER(iwp) ::  k !<
+
+       REAL(wp)    ::  dqdt_precip !<
+
+       DO  k = nzb_2d(j,i)+1, nzt
+
+          IF ( qc_1d(k) > ql_crit )  THEN
+             dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
+          ELSE
+             dqdt_precip = 0.0_wp
+          ENDIF
+
+          qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro
+          q_1d(k)  = q_1d(k)  - dqdt_precip * dt_micro
+          pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k)
+
+!
+!--       Compute the rain rate
+          prr(nzb_2d(j,i)+1,j,i) = prr(nzb_2d(j,i)+1,j,i) +                    &
+                                   dqdt_precip * dzw(k)
+
+       ENDDO
+
+    END SUBROUTINE autoconversion_kessler_ij
 
 !------------------------------------------------------------------------------!
@@ -1384 +1499 @@
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nzb_s_inner, nzt
 
        USE kinds
@@ -1398 +1513 @@
        REAL(wp)     ::  phi_ac            !<
        REAL(wp)     ::  tau_cloud         !<
-       REAL(wp)     ::  xc                !<
 
        DO  k = nzb_s_inner(j,i)+1, nzt
@@ -1453 +1567 @@
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nzb_s_inner, nzt
 
        USE kinds
@@ -1507 +1621 @@
        USE cloud_parameters,                                                   &
            ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
-                  dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,   &
-                  l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,      &
+                  dpirho_l, eps_sb, hyrho, kin_vis_air, l_d_cp, l_d_r,         &
+                  l_v, r_v, schmidt_p_1d3, thermal_conductivity_l,             &
                   t_d_pt, ventilation_effect
 
@@ -1518 +1632 @@
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nzb_s_inner, nzt
 
        USE kinds
@@ -1655 +1769 @@
            ONLY:  eps_sb, hyrho, l_d_cp, prr, pt_d_t, sed_qc_const
 
-       USE constants,                                                          &
-           ONLY:  pi
-
        USE control_parameters,                                                 &
-           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
-                  intermediate_timestep_count, precipitation
+           ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
+                  intermediate_timestep_count
 
        USE indices,                                                            &
@@ -1701 +1812 @@
 !
 !--       Compute the precipitation rate of cloud (fog) droplets
-          IF ( precipitation )  THEN
-             IF ( call_microphysics_at_all_substeps )  THEN
-                prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *               &
+          IF ( call_microphysics_at_all_substeps )  THEN
+             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *               &
                              weight_substep(intermediate_timestep_count)
-             ELSE
-                prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)
-             ENDIF
+          ELSE
+             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
           ENDIF
 
@@ -1727 +1836 @@
 
        USE cloud_parameters,                                                   &
-           ONLY:  a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,        &
-                  limiter_sedimentation, l_d_cp, precipitation_amount, prr,    &
-                  pt_d_t, stp
+           ONLY:  a_term, b_term, c_term, dpirho_l, eps_sb, hyrho,             &
+                  limiter_sedimentation, l_d_cp, prr, pt_d_t
 
        USE control_parameters,                                                 &
-           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
-                  dt_3d, intermediate_timestep_count,                          &
-                  intermediate_timestep_count_max,                             &
-                  precipitation_amount_interval, time_do2d_xy
+           ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
+                  intermediate_timestep_count
 
        USE indices,                                                            &
@@ -1757 +1863 @@
        REAL(wp)     ::  d_min                      !<
        REAL(wp)     ::  dr                         !<
-       REAL(wp)     ::  dt_sedi                    !<
        REAL(wp)     ::  flux                       !<
        REAL(wp)     ::  lambda_r                   !<
@@ -1765 +1870 @@
        REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
        REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr     !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_qr     !<
        REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
        REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
@@ -1774 +1877 @@
        REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
 
-       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
 !
 !--    Compute velocities 
@@ -1939 +2041 @@
            ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
                   intermediate_timestep_count, intermediate_timestep_count_max,&
-                  precipitation_amount_interval,         &
-                  time_do2d_xy
+                  precipitation_amount_interval, time_do2d_xy
 
        USE indices,                                                            &

palm/trunk/SOURCE/mod_particle_attributes.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
+! +collision_algorithm, all_or_nothing, average_impact
 ! 
-! 
+! Tails removed.
+!
 ! Former revisions:
 ! -----------------
@@ -47 +49 @@
     USE kinds
 
-    CHARACTER(LEN=15) ::  bc_par_lr = 'cyclic',  bc_par_ns = 'cyclic',         &
-                          bc_par_b  = 'reflect', bc_par_t  = 'absorb',         &
-                          collision_kernel = 'none'
+    CHARACTER(LEN=15) ::  bc_par_lr = 'cyclic'                    !< left/right boundary condition
+    CHARACTER(LEN=15) ::  bc_par_ns = 'cyclic'                    !< north/south boundary condition
+    CHARACTER(LEN=15) ::  bc_par_b  = 'reflect'                   !< bottom boundary condition
+    CHARACTER(LEN=15) ::  bc_par_t  = 'absorb'                    !< top boundary condition
+    CHARACTER(LEN=15) ::  collision_algorithm = 'all_or_nothing'  !< collision algorithm
+    CHARACTER(LEN=15) ::  collision_kernel = 'none'               !< collision kernel
 
-    INTEGER(iwp) ::  deleted_particles = 0, deleted_tails = 0,                 &
+    INTEGER(iwp) ::  deleted_particles = 0,                                    &
                      dissipation_classes = 10, ibc_par_lr,                     &
                      ibc_par_ns, ibc_par_b, ibc_par_t, iran_part = -1234567,   &
                      maximum_number_of_particles = 0,                          &
-                     maximum_number_of_tailpoints = 100,                       &
-                     maximum_number_of_tails = 0,                              &
                      min_nr_particle = 50,                                     &
                      mpi_particle_type,                                        &
                      number_of_particles = 0,                                  &
-                     number_of_particle_groups = 1, number_of_tails = 0,       &
-                     number_of_initial_tails = 0, number_of_sublayers = 20,    &
+                     number_of_particle_groups = 1,                            &
+                     number_of_sublayers = 20,                                 &
                      offset_ocean_nzt = 0,                                     &
                      offset_ocean_nzt_m1 = 0, particles_per_point = 1,         &
                      particle_file_count = 0, radius_classes = 20,             &
-                     skip_particles_for_tail = 100, sort_count = 0,            &
-                     total_number_of_particles, total_number_of_tails = 0,     &
+                     sort_count = 0,                                           &
+                     total_number_of_particles,                                &
                      trlp_count_sum, trlp_count_recv_sum, trrp_count_sum,      &
                      trrp_count_recv_sum, trsp_count_sum, trsp_count_recv_sum, &
@@ -73 +76 @@
     INTEGER(iwp), PARAMETER ::  max_number_of_particle_groups = 10
 
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE     ::  new_tail_id
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  prt_count
 
-    LOGICAL ::  hall_kernel = .FALSE., palm_kernel = .FALSE.,                  &
+    LOGICAL ::  all_or_nothing = .FALSE., average_impact = .FALSE.,            &
+                hall_kernel = .FALSE., palm_kernel = .FALSE.,                  &
                 particle_advection = .FALSE., random_start_position = .FALSE., &
                 read_particles_from_restartfile = .TRUE.,                      &
                 seed_follows_topography = .FALSE.,                             &
                 uniform_particles = .TRUE., use_kernel_tables = .FALSE.,       &
-                use_particle_tails = .FALSE., use_sgs_for_particles = .FALSE., &
-                wang_kernel = .FALSE., write_particle_statistics = .FALSE.
+                use_sgs_for_particles = .FALSE., wang_kernel = .FALSE.,        &
+                write_particle_statistics = .FALSE.
 
     LOGICAL, DIMENSION(max_number_of_particle_groups) ::                       &
                 vertical_particle_advection = .TRUE.
 
-    LOGICAL, DIMENSION(:), ALLOCATABLE ::  tail_mask
-
     REAL(wp) ::  alloc_factor = 20.0_wp, c_0 = 3.0_wp,                         &
                  dt_min_part = 0.0002_wp, dt_prel = 9999999.9_wp,              &
                  dt_write_particle_data = 9999999.9_wp,                        &
-                 dvrp_psize = 9999999.9_wp, end_time_prel = 9999999.9_wp,      &
+                 end_time_prel = 9999999.9_wp,                                 &
                  initial_weighting_factor = 1.0_wp,                            &
-                 maximum_tailpoint_age = 100000.0_wp,                          &
-                 minimum_tailpoint_distance = 0.0_wp,                          &
                  particle_advection_start = 0.0_wp,                            &
                  sgs_wfu_part = 0.3333333_wp, sgs_wfv_part = 0.3333333_wp,     &
@@ -108 +107 @@
                  pss = 9999999.9_wp, pst = 9999999.9_wp, radius = 9999999.9_wp
 
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  particle_tail_coordinates
-
     REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0
 
     TYPE particle_type
         SEQUENCE
-#if defined( __twocachelines )
-        REAL(wp)     ::  radius
-        REAL(sp)     ::  x, y, z, speed_x, speed_y, speed_z
-        REAL(wp)     ::  weight_factor, rvar1, dt_sum
-        INTEGER(iwp) ::  class
-        LOGICAL      ::  particle_mask
-
-        REAL(wp)     ::  dvrp_psize, rvar2, rvar3
-        REAL(sp)     ::  age, origin_x, origin_y, origin_z, e_m, age_m
-        INTEGER(iwp) ::  group, tailpoints, tail_id, block_nr
-#else
         REAL(wp)     ::  radius, age, age_m, dt_sum, dvrp_psize, e_m,          &
                          origin_x, origin_y, origin_z, rvar1, rvar2, rvar3,    &
@@ -131 +117 @@
         LOGICAL      ::  particle_mask
         INTEGER(iwp) ::  block_nr
-#endif
-!-- One 32 Bit word for 64 Bit alignment in Type declaration
     END TYPE particle_type
 
-    TYPE(particle_type), DIMENSION(:), POINTER             ::  particles
-    TYPE(particle_type)                                    ::  zero_particle
+    TYPE(particle_type), DIMENSION(:), POINTER ::  particles
+    TYPE(particle_type)                        ::  zero_particle
 
     TYPE particle_groups_type

palm/trunk/SOURCE/modules.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-!
+! icloud_scheme removed. microphysics_sat_adjust, microphysics_kessler, 
+! microphysics_seifert added. 
 !
 ! Former revisions:
@@ -626 +627 @@
                              bc_sa_t = 'neumann', &
                              bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', &
-                             cloud_scheme = 'seifert_beheng', &
+                             cloud_scheme = 'saturation_adjust', &
                              coupling_mode = 'uncoupled', &
                              coupling_mode_remote = 'uncoupled', &
@@ -665 +666 @@
                      grid_level, ibc_e_b, ibc_p_b, ibc_p_t, &
                      ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, &
-                     ibc_sa_t, ibc_uv_b, ibc_uv_t, icloud_scheme, &
+                     ibc_sa_t, ibc_uv_b, ibc_uv_t, &
                      inflow_disturbance_begin = -1, inflow_disturbance_end = -1, &
                      intermediate_timestep_count, intermediate_timestep_count_max, &
@@ -734 +735 @@
                 lsf_vert = .TRUE., lptnudge = .FALSE., lqnudge = .FALSE., &
                 lunudge = .FALSE., lvnudge = .FALSE., lwnudge = .FALSE., &
-                masking_method = .FALSE., mg_switch_to_pe0 = .FALSE., &
+                masking_method = .FALSE., &
+                microphysics_sat_adjust = .FALSE., &
+                microphysics_kessler = .FALSE., &
+                microphysics_seifert = .FALSE., &
+                mg_switch_to_pe0 = .FALSE., &
                 monotonic_adjustment = .FALSE.
     LOGICAL ::  nest_bound_l = .FALSE. !< nested boundary on left side

palm/trunk/SOURCE/package_parin.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! +collision_algorithm
+!
+! Tails removed.
 ! 
 ! Former revisions:
@@ -137 +139 @@
     USE particle_attributes,                                                   &
         ONLY:  alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,         &
-               collision_kernel, density_ratio, dissipation_classes,           &
-               dt_min_part, dt_prel, dt_write_particle_data, dvrp_psize,       &
+               collision_algorithm, collision_kernel, density_ratio,           &
+               dissipation_classes, dt_min_part, dt_prel,                      &
+               dt_write_particle_data,                                         &
                end_time_prel, initial_weighting_factor,                        &
-               maximum_number_of_tailpoints,                                   &
-               maximum_tailpoint_age, minimum_tailpoint_distance,              &
                min_nr_particle, number_of_particle_groups, particles_per_point,&
                particle_advection, particle_advection_start, pdx, pdy, pdz,    &
                psb, psl, psn, psr, pss, pst, radius, radius_classes,           &
                random_start_position, read_particles_from_restartfile,         &
-               seed_follows_topography, skip_particles_for_tail,               &
-               use_particle_tails, use_sgs_for_particles,                      &
+               seed_follows_topography, use_sgs_for_particles,                 &
                vertical_particle_advection, write_particle_statistics
 
@@ -186 +186 @@
                                   dt_dvrp, dvrpsize_interval, dvrp_directory,  &
                                   dvrp_file, dvrp_host, dvrp_output,           &
-                                  dvrp_password, dvrp_psize, dvrp_username,    &
+                                  dvrp_password, dvrp_username,                &
                                   groundplate_color, isosurface_color,         &
                                   mode_dvrp, particle_color, particle_dvrpsize,&
@@ -199 +199 @@
 
     NAMELIST /particles_par/      alloc_factor, bc_par_b, bc_par_lr,           &
-                                  bc_par_ns, bc_par_t, collision_kernel,       &
+                                  bc_par_ns, bc_par_t, collision_algorithm,    &
+                                  collision_kernel,                            &
                                   density_ratio, dissipation_classes, dt_dopts,&
                                   dt_min_part, dt_prel,                        &
                                   dt_write_particle_data,                      &
                                   end_time_prel, initial_weighting_factor,     &
-                                  maximum_number_of_tailpoints,                &
-                                  maximum_tailpoint_age,                       &
-                                  minimum_tailpoint_distance,                  &
                                   min_nr_particle,                             &
                                   number_of_particle_groups,                   &
@@ -216 +214 @@
                                   read_particles_from_restartfile,             &
                                   seed_follows_topography,                     &
-                                  skip_particles_for_tail,                     &
-                                  use_particle_tails, use_sgs_for_particles,   &
+                                  use_sgs_for_particles,                       &
                                   vertical_particle_advection,                 &
                                   write_particle_statistics

palm/trunk/SOURCE/prognostic_equations.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Kessler microphysics scheme moved to microphysics.  
 ! 
 ! Former revisions:
@@ -190 +190 @@
     USE control_parameters,                                                    &
         ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
-               cloud_top_radiation, constant_diffusion, dp_external,       &
+               cloud_top_radiation, constant_diffusion, dp_external,           &
                dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
-               icloud_scheme, inflow_l, intermediate_timestep_count,           &
+               inflow_l, intermediate_timestep_count,                          &
                intermediate_timestep_count_max, large_scale_forcing,           &
-               large_scale_subsidence, neutral, nudging, ocean, outflow_l,     &
-               outflow_s, passive_scalar, precipitation,                       &
-               prho_reference, prho_reference, prho_reference, pt_reference,   &
-               pt_reference, pt_reference, scalar_advec, scalar_advec,         &
-               simulated_time, sloping_surface,                                &
+               large_scale_subsidence, microphysics_seifert,                   &
+               microphysics_sat_adjust, neutral, nudging, ocean, outflow_l,    &
+               outflow_s, passive_scalar, prho_reference, prho_reference,      &
+               prho_reference, pt_reference, pt_reference, pt_reference,       &
+               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
                timestep_scheme, tsc, use_subsidence_tendencies,                &
                use_upstream_for_tke, wall_heatflux,                            &
@@ -247 +247 @@
     USE buoyancy_mod,                                                          &
         ONLY:  buoyancy, buoyancy_acc
-
-    USE calc_precipitation_mod,                                                &
-        ONLY:  calc_precipitation
 
     USE calc_radiation_mod,                                                    &
@@ -271 +268 @@
     USE diffusion_w_mod,                                                       &
         ONLY:  diffusion_w, diffusion_w_acc
-
-    USE impact_of_latent_heat_mod,                                             &
-        ONLY:  impact_of_latent_heat
 
     USE kinds
@@ -377 +371 @@
        DO  j = nys, nyn
 !
-!--       If required, calculate cloud microphysical impacts (two-moment scheme)
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
+!--       If required, calculate cloud microphysics
+          IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.      &
                ( intermediate_timestep_count == 1  .OR.                        &
                  call_microphysics_at_all_substeps )                           &
@@ -584 +578 @@
              IF ( cloud_top_radiation )  THEN
                 CALL calc_radiation( i, j )
-             ENDIF
-
-!
-!--          If required compute impact of latent heat due to precipitation
-             IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
-                  precipitation )  THEN
-                CALL impact_of_latent_heat( i, j ) 
              ENDIF
 
@@ -729 +716 @@
 
 !
-!--          If required compute decrease of total water content due to
-!--          precipitation
-             IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
-                  precipitation )  THEN
-                CALL calc_precipitation( i, j )
-             ENDIF
-
-!
 !--          Sink or source of scalar concentration due to canopy elements
              IF ( plant_canopy )  CALL plant_canopy_model( i, j, 5 )
@@ -788 +767 @@
 !--          If required, calculate prognostic equations for rain water content 
 !--          and rain drop concentration
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
-                  precipitation )  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 !
 !--             Calculate prognostic equation for rain water content
@@ -807 +785 @@
                 CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
 
-                CALL user_actions( i, j, 'qr-tendency' )
 !
 !--             Prognostic equation for rain water content
@@ -848 +825 @@
                 CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux )
 
-                CALL user_actions( i, j, 'nr-tendency' )
 !
 !--             Prognostic equation for rain drop concentration
@@ -971 +947 @@
 
 !
-!-- If required, calculate cloud microphysical impacts (two-moment scheme)
-    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                       &
+!-- If required, calculate cloud microphysical impacts
+    IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.            &
          ( intermediate_timestep_count == 1  .OR.                              &
            call_microphysics_at_all_substeps )                                 &
@@ -1254 +1230 @@
 
 !
-!--    If required compute impact of latent heat due to precipitation
-       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
-          CALL impact_of_latent_heat
-       ENDIF
-
-!
 !--    Consideration of heat sources within the plant canopy
        IF ( plant_canopy .AND. ( cthf /= 0.0_wp ) ) THEN
@@ -1459 +1429 @@
        
 !
-!--    If required compute decrease of total water content due to
-!--    precipitation
-       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
-          CALL calc_precipitation
-       ENDIF
-
-!
 !--    Sink or source of scalar concentration due to canopy elements
        IF ( plant_canopy ) CALL plant_canopy_model( 5 )
@@ -1530 +1493 @@
 !--    If required, calculate prognostic equations for rain water content 
 !--    and rain drop concentration
-       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 
           CALL cpu_log( log_point(52), 'qr-equation', 'start' )
@@ -1565 +1528 @@
 
           CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
-
-          CALL user_actions( 'qr-tendency' )
 
 !
@@ -1639 +1600 @@
 
           CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
-
-          CALL user_actions( 'nr-tendency' )
 
 !
@@ -1821 +1780 @@
 !
 !-- If required, calculate cloud microphysical impacts (two-moment scheme)
-    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                       &
+    IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.            &
          ( intermediate_timestep_count == 1  .OR.                              &
            call_microphysics_at_all_substeps )                                 &
@@ -2079 +2038 @@
 
 !
-!--    If required compute impact of latent heat due to precipitation
-       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
-          CALL impact_of_latent_heat
-       ENDIF
-
-!
 !--    Consideration of heat sources within the plant canopy
        IF ( plant_canopy .AND. ( cthf /= 0.0_wp ) ) THEN
@@ -2256 +2209 @@
 
 !
-!--    If required compute decrease of total water content due to
-!--    precipitation
-       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
-          CALL calc_precipitation
-       ENDIF
-
-!
 !--    Sink or source of scalar concentration due to canopy elements
        IF ( plant_canopy ) CALL plant_canopy_model( 5 )
@@ -2311 +2257 @@
 !--    If required, calculate prognostic equations for rain water content 
 !--    and rain drop concentration
-       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 
           CALL cpu_log( log_point(52), 'qr-equation', 'start' )
@@ -2346 +2292 @@
           CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
 
-          CALL user_actions( 'qr-tendency' )
-
 !
 !--       Prognostic equation for rain water content
@@ -2403 +2347 @@
 
           CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
-
-          CALL user_actions( 'nr-tendency' )
 
 !

palm/trunk/SOURCE/surface_layer_fluxes.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! icloud_scheme replaced by microphysics_*
 ! 
 ! Former revisions:
@@ -134 +134 @@
     USE control_parameters,                                                    &
         ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
-               coupling_mode, g, humidity, ibc_e_b, ibc_pt_b, icloud_scheme,   &
+               coupling_mode, g, humidity, ibc_e_b, ibc_pt_b,                  &
                initializing_actions, kappa, intermediate_timestep_count,       &
                intermediate_timestep_count_max, large_scale_forcing, lsf_surf, &
-               message_string, most_method, neutral, passive_scalar,           &
-               precipitation, pt_surface, q_surface, run_coupled,              &
+               message_string, microphysics_seifert, most_method, neutral,     &
+               passive_scalar, pt_surface, q_surface, run_coupled,             &
                surface_pressure, simulated_time, terminate_run, zeta_max,      &
                zeta_min 
@@ -906 +906 @@
 !
 !--    If required compute qr* and nr*
-       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )   &
+       IF ( cloud_physics  .AND.  microphysics_seifert )   &
        THEN
 
@@ -1040 +1040 @@
 !
 !--    Compute (turbulent) fluxes of rain water content and rain drop conc.
-       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                    &
-            precipitation )  THEN
+       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
           !$OMP PARALLEL DO
           !$acc kernels loop independent present( nrs, nrsws, qrs, qrsws, us )

palm/trunk/SOURCE/swap_timelevel.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-!
+! icloud_scheme replaced by microphysics_*
 !
 ! Former revisions:
@@ -97 +97 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_physics, constant_diffusion, humidity, icloud_scheme,     &
-               neutral, ocean, passive_scalar, precipitation, timestep_count
+        ONLY:  cloud_physics, constant_diffusion, humidity,                    &
+               microphysics_seifert, neutral, ocean, passive_scalar,           &
+               timestep_count
 
     USE indices,                                                               &
@@ -160 +161 @@
     IF ( humidity  .OR.  passive_scalar )  THEN
        q = q_p             
-       IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
           qr = qr_p
           nr = nr_p
@@ -193 +194 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              q => q_1;    q_p => q_2
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
-                  precipitation )  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 qr => qr_1;    qr_p => qr_2
                 nr => nr_1;    nr_p => nr_2
@@ -222 +222 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              q => q_2;    q_p => q_1
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
-                  precipitation)  THEN
+             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 qr => qr_2;    qr_p => qr_1
                 nr => nr_2;    nr_p => nr_1

palm/trunk/SOURCE/time_integration.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-!
+! icloud_scheme replaced by microphysics_*
 !
 ! Former revisions:
@@ -213 +213 @@
                dt_dopr_listing, dt_dosp, dt_dots, dt_dvrp, dt_run_control,     &
                end_time, first_call_lpm, galilei_transformation, humidity,     &
-               icloud_scheme, intermediate_timestep_count,                     &
+               intermediate_timestep_count,                                    &
                intermediate_timestep_count_max, large_scale_forcing,           &
-               loop_optimization, lsf_surf, lsf_vert, masks, mid, nest_domain, &
+               loop_optimization, lsf_surf, lsf_vert, masks,                   &
+               microphysics_seifert, mid, nest_domain,                         &
                neutral, nr_timesteps_this_run, nudging,                        &
-               ocean, on_device, passive_scalar, precipitation,                &
+               ocean, on_device, passive_scalar,                               &
                prho_reference, pt_reference, pt_slope_offset, random_heatflux, &
                run_coupled, simulated_time, simulated_time_chr,                &
@@ -486 +487 @@
                 IF (humidity  .OR.  passive_scalar)  THEN
                    CALL exchange_horiz( q_p, nbgp )
-                   IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                   IF ( cloud_physics .AND. microphysics_seifert )  THEN
                       CALL exchange_horiz( qr_p, nbgp )
                       CALL exchange_horiz( nr_p, nbgp )
@@ -553 +554 @@
                 IF (humidity  .OR.  passive_scalar)  THEN
                    CALL exchange_horiz( q_p, nbgp )
-                   IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                   IF ( cloud_physics .AND. microphysics_seifert )  THEN
                       CALL exchange_horiz( qr_p, nbgp )
                       CALL exchange_horiz( nr_p, nbgp )
@@ -642 +643 @@
              IF (humidity  .OR.  passive_scalar)  THEN
                 CALL exchange_horiz( q_p, nbgp )
-                IF ( cloud_physics .AND. icloud_scheme == 0  .AND.  &
-                     precipitation )  THEN
+                IF ( cloud_physics .AND. microphysics_seifert )  THEN
                    CALL exchange_horiz( qr_p, nbgp )
                    CALL exchange_horiz( nr_p, nbgp )

palm/trunk/SOURCE/user_actions.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! qr/nr-tendency removed.
 ! 
 ! Former revisions:
@@ -161 +161 @@
 
 
-          CASE ( 'qr-tendency' )
-
-
-          CASE ( 'nr-tendency' )
-
-
           CASE DEFAULT
              message_string = 'unknown location "' // location // '"'
@@ -225 +219 @@
 
 
-          CASE ( 'qr-tendency' )
-
-
-          CASE ( 'nr-tendency' )
-
-
           CASE ( 'before_timestep', 'after_integration', 'after_timestep' )
              message_string = 'location "' // location // '" is not ' // &

palm/trunk/SOURCE/write_3d_binary.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! icloud_scheme replaced by microphysics_*
 ! 
 ! Former revisions:
@@ -102 +102 @@
     USE control_parameters,                                                    &
         ONLY:  iran, humidity, passive_scalar, cloud_physics, cloud_droplets,  &
-               icloud_scheme, precipitation, ocean, topography
+               microphysics_seifert, ocean, topography
                 
     USE indices,                                                               &
@@ -240 +240 @@
              WRITE ( 14 )  'ql_av               ';  WRITE ( 14 )  ql_av
           ENDIF
-          IF ( icloud_scheme == 0  .AND.  .NOT. cloud_droplets )  THEN
+          IF ( .NOT. cloud_droplets )  THEN
              WRITE ( 14 )  'qc                  ';  WRITE ( 14 )  qc
              IF ( ALLOCATED( qc_av ) )  THEN
                 WRITE ( 14 )  'qc_av               ';  WRITE ( 14 )  qc_av
              ENDIF
-             IF ( precipitation )  THEN
+             IF ( microphysics_seifert )  THEN
                 WRITE ( 14 )  'nr                  ';  WRITE ( 14 )  nr
                 IF ( ALLOCATED( nr_av ) )  THEN