Changeset 849 for palm/trunk


Ignore:
Timestamp:
Mar 15, 2012 10:35:09 AM (10 years ago)
Author:
raasch
Message:

Changed:


Original routine advec_particles split into several new subroutines and renamed
lpm.
init_particles renamed lpm_init
user_advec_particles renamed user_lpm_advec,
particle_boundary_conds renamed lpm_boundary_conds,
set_particle_attributes renamed lpm_set_attributes,
user_init_particles renamed user_lpm_init,
user_particle_attributes renamed user_lpm_set_attributes
(Makefile, lpm_droplet_collision, lpm_droplet_condensation, init_3d_model, modules, palm, read_var_list, time_integration, write_var_list, deleted: advec_particles, init_particles, particle_boundary_conds, set_particle_attributes, user_advec_particles, user_init_particles, user_particle_attributes, new: lpm, lpm_advec, lpm_boundary_conds, lpm_calc_liquid_water_content, lpm_data_output_particles, lpm_droplet_collision, lpm_drollet_condensation, lpm_exchange_horiz, lpm_extend_particle_array, lpm_extend_tails, lpm_extend_tail_array, lpm_init, lpm_init_sgs_tke, lpm_pack_arrays, lpm_read_restart_file, lpm_release_set, lpm_set_attributes, lpm_sort_arrays, lpm_write_exchange_statistics, lpm_write_restart_file, user_lpm_advec, user_lpm_init, user_lpm_set_attributes

Location:
palm/trunk/SOURCE
Files:
16 added
9 edited
7 moved

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r829 r849  
    44# Current revisions:
    55# ------------------
    6 #
     6# +lpm_advec, lpm_data_output_particles, lpm_droplet_collision,
     7# lpm_droplet_condensation,
     8# lpm_exchange_horiz, lpm_extend_particle_array, lpm_extend_tails,
     9# lpm_extend_tail_array, lpm_init_sgs_tke, lpm_pack_arrays,
     10# lpm_read_restart_file, lpm_release_set, lpm_sort_arrays,
     11# lpm_write_exchange_statistics, lpm_write_restart_file
     12# advec_particles renamed lpm,
     13# init_particles renamed lpm_init,
     14# user_advec_particles renamed user_lpm_advec,
     15# particle_boundary_conds renamed lpm_boundary_conds,
     16# set_particle_attributes renamed lpm_set_attributes
     17# user_init_particles renamed user_lpm_init,
     18# user_particle_attributes renamed user_lpm_set_attributes
    719#
    820# Former revisions:
     
    6880PROG =  palm
    6981
    70 RCS = advec_particles.f90 advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \
     82RCS =   advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \
    7183        advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \
    7284        advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \
     
    8799        inflow_turbulence.f90 init_1d_model.f90 init_3d_model.f90 \
    88100        init_advec.f90 init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 \
    89         init_grid.f90 init_masks.f90 init_ocean.f90 init_particles.f90 \
     101        init_grid.f90 init_masks.f90 init_ocean.f90 \
    90102        init_pegrid.f90 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
    91103        interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \
    92104        local_stop.f90 local_system.f90 local_tremain.f90 \
    93         local_tremain_ini.f90 lpm_collision_kernels.f90 message.f90 \
    94         modules.f90 netcdf.f90 package_parin.f90 palm.f90 parin.f90 \
    95         particle_boundary_conds.f90 plant_canopy_model.f90 poisfft.f90 \
     105        local_tremain_ini.f90 lpm.f90 lpm_advec.f90 lpm_boundary_conds.f90 \
     106        lpm_calc_liquid_water_content.f90 lpm_collision_kernels.f90 \
     107        lpm_data_output_particles.f90 lpm_droplet_collision.f90 \
     108        lpm_droplet_condensation.f90 lpm_exchange_horiz.f90 \
     109        lpm_extend_particle_array.f90 lpm_extend_tails.f90 \
     110        lpm_extend_tail_array.f90 lpm_init.f90 lpm_init_sgs_tke.f90 \
     111        lpm_pack_arrays.f90 lpm_read_restart_file.f90 lpm_release_set.f90 \
     112        lpm_set_attributes.f90 lpm_sort_arrays.f90 \
     113        lpm_write_exchange_statistics.f90 lpm_write_restart_file.f90 \
     114        message.f90 modules.f90 netcdf.f90 package_parin.f90 palm.f90 \
     115        parin.f90 plant_canopy_model.f90 poisfft.f90 \
    96116        poisfft_hybrid.f90 poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \
    97117        production_e.f90 prognostic_equations.f90 random_function.f90 \
    98118        random_gauss.f90 read_3d_binary.f90 read_var_list.f90 run_control.f90 \
    99         set_particle_attributes.f90 set_slicer_attributes_dvrp.f90 \
    100         singleton.f90 sor.f90 spline_x.f90 \
    101         spline_y.f90 spline_z.f90 subsidence.f90 \
    102         sum_up_3d_data.f90 surface_coupler.f90 \
    103         swap_timelevel.f90 temperton_fft.f90 time_integration.f90 \
    104         time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \
    105         transpose.f90 user_3d_data_averaging.f90 user_actions.f90 \
    106         user_additional_routines.f90 user_advec_particles.f90 \
     119        set_slicer_attributes_dvrp.f90 singleton.f90 sor.f90 spline_x.f90 \
     120        spline_y.f90 spline_z.f90 subsidence.f90 sum_up_3d_data.f90 \
     121        surface_coupler.f90 swap_timelevel.f90 temperton_fft.f90 \
     122        time_integration.f90 time_to_string.f90 timestep.f90 \
     123        timestep_scheme_steering.f90 transpose.f90 user_3d_data_averaging.f90 \
     124        user_actions.f90 user_additional_routines.f90 \
    107125        user_check_data_output.f90 user_check_data_output_pr.f90 \
    108126        user_check_parameters.f90 user_data_output_2d.f90 \
     
    110128        user_data_output_mask.f90 user_define_netcdf_grid.f90 \
    111129        user_dvrp_coltab.f90 user_header.f90 user_init.f90 \
    112         user_init_3d_model.f90 user_init_grid.f90 user_init_particles.f90 \
    113         user_init_plant_canopy.f90 user_last_actions.f90 user_module.f90 \
    114         user_parin.f90 user_particle_attributes.f90 user_read_restart_data.f90 \
     130        user_init_3d_model.f90 user_init_grid.f90 \
     131        user_init_plant_canopy.f90 user_last_actions.f90 user_lpm_advec.f90 \
     132        user_lpm_init.f90 user_lpm_set_attributes.f90 user_module.f90 \
     133        user_parin.f90 user_read_restart_data.f90 \
    115134        user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
    116135        write_3d_binary.f90 write_compressed.f90 write_var_list.f90
    117136
    118 OBJS =  advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \
    119         advec_s_ups.o advec_u_pw.o advec_u_up.o advec_u_ups.o \
    120         advec_ws.o advec_v_pw.o advec_v_up.o advec_v_ups.o advec_w_pw.o \
    121         advec_w_up.o advec_w_ups.o asselin_filter.o average_3d_data.o \
    122         boundary_conds.o buoyancy.o calc_liquid_water_content.o \
    123         calc_precipitation.o calc_radiation.o calc_spectra.o \
    124         check_for_restart.o check_open.o check_parameters.o close_file.o \
    125         compute_vpt.o coriolis.o cpu_log.o cpu_statistics.o data_log.o \
    126         data_output_dvrp.o data_output_mask.o data_output_profiles.o \
    127         data_output_ptseries.o \
     137OBJS =  advec_s_bc.o advec_s_pw.o advec_s_up.o advec_s_ups.o advec_u_pw.o \
     138        advec_u_up.o advec_u_ups.o advec_ws.o advec_v_pw.o advec_v_up.o \
     139        advec_v_ups.o advec_w_pw.o advec_w_up.o advec_w_ups.o asselin_filter.o \
     140        average_3d_data.o boundary_conds.o buoyancy.o \
     141        calc_liquid_water_content.o calc_precipitation.o calc_radiation.o \
     142        calc_spectra.o check_for_restart.o check_open.o check_parameters.o \
     143        close_file.o compute_vpt.o coriolis.o cpu_log.o cpu_statistics.o \
     144        data_log.o data_output_dvrp.o data_output_mask.o \
     145        data_output_profiles.o data_output_ptseries.o \
    128146        data_output_spectra.o data_output_tseries.o data_output_2d.o \
    129147        data_output_3d.o diffusion_e.o diffusion_s.o diffusion_u.o \
     
    133151        header.o impact_of_latent_heat.o inflow_turbulence.o init_1d_model.o \
    134152        init_3d_model.o init_advec.o init_cloud_physics.o init_coupling.o \
    135         init_dvrp.o init_grid.o init_masks.o init_ocean.o init_particles.o \
    136         init_pegrid.o init_pt_anomaly.o init_rankine.o init_slope.o \
     153        init_dvrp.o init_grid.o init_masks.o init_ocean.o init_pegrid.o \
     154        init_pt_anomaly.o init_rankine.o init_slope.o \
    137155        interaction_droplets_ptq.o local_flush.o local_getenv.o local_stop.o \
    138         local_system.o local_tremain.o local_tremain_ini.o \
    139         lpm_collision_kernels.f90 message.o modules.o \
    140         netcdf.o package_parin.o palm.o parin.o particle_boundary_conds.o \
    141         plant_canopy_model.o poisfft.o \
     156        local_system.o local_tremain.o local_tremain_ini.o lpm.o lpm_advec.o \
     157        lpm_boundary_conds.o lpm_calc_liquid_water_content.o \
     158        lpm_collision_kernels.o lpm_data_output_particles.o \
     159        lpm_droplet_collision.o lpm_droplet_condensation.o \
     160        lpm_exchange_horiz.o lpm_extend_particle_array.o lpm_extend_tails.o \
     161        lpm_extend_tail_array.o lpm_init.o lpm_init_sgs_tke.o \
     162        lpm_pack_arrays.o lpm_read_restart_file.o lpm_release_set.o \
     163        lpm_set_attributes.o lpm_sort_arrays.o lpm_write_exchange_statistics.o \
     164        lpm_write_restart_file.o message.o modules.o netcdf.o package_parin.o \
     165        palm.o parin.o plant_canopy_model.o poisfft.o \
    142166        poisfft_hybrid.o poismg.o prandtl_fluxes.o pres.o print_1d.o \
    143167        production_e.o prognostic_equations.o random_function.o random_gauss.o \
    144168        read_3d_binary.o read_var_list.o run_control.o \
    145         set_particle_attributes.o set_slicer_attributes_dvrp.o \
    146         singleton.o sor.o spline_x.o spline_y.o \
    147         spline_z.o subsidence.o sum_up_3d_data.o surface_coupler.o \
    148         swap_timelevel.o \
    149         temperton_fft.o time_integration.o time_to_string.o \
    150         timestep.o timestep_scheme_steering.o transpose.o \
     169        set_slicer_attributes_dvrp.o singleton.o sor.o spline_x.o spline_y.o \
     170        spline_z.o subsidence.o sum_up_3d_data.o surface_coupler.o \
     171        swap_timelevel.o temperton_fft.o time_integration.o time_to_string.o \
     172        timestep.o timestep_scheme_steering.o transpose.o \
    151173        user_3d_data_averaging.o user_actions.o user_additional_routines.o \
    152         user_advec_particles.o user_check_data_output.o \
    153         user_check_data_output_pr.o user_check_parameters.o \
    154         user_data_output_2d.o user_data_output_3d.o user_data_output_mask.o user_data_output_dvrp.o \
     174        user_check_data_output.o user_check_data_output_pr.o \
     175        user_check_parameters.o user_data_output_2d.o user_data_output_3d.o \
     176        user_data_output_mask.o user_data_output_dvrp.o \
    155177        user_define_netcdf_grid.o user_dvrp_coltab.o user_header.o \
    156178        user_init.o user_init_3d_model.o user_init_grid.o \
    157         user_init_particles.o user_init_plant_canopy.o user_last_actions.o \
    158         user_module.o user_parin.o user_particle_attributes.o \
     179        user_init_plant_canopy.o user_last_actions.o user_lpm_advec.o \
     180        user_lpm_init.o user_lpm_set_attributes.o user_module.o user_parin.o \
    159181        user_read_restart_data.o user_spectra.o user_statistics.o \
    160         wall_fluxes.o write_3d_binary.o write_compressed.o \
    161         write_var_list.o
     182        wall_fluxes.o write_3d_binary.o write_compressed.o write_var_list.o
    162183
    163184CC = cc
     
    188209
    189210
    190 advec_particles.o: modules.o lpm_collision_kernels.o random_function.o
     211advec_particles.o: modules.o
    191212advec_s_bc.o: modules.o
    192213advec_s_pw.o: modules.o
     
    254275init_masks.o: modules.o
    255276init_ocean.o: modules.o eqn_state_seawater.o
    256 init_particles.o: modules.o lpm_collision_kernels.o random_function.o
    257277init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o
    258278init_pt_anomaly.o: modules.o
     
    264284local_tremain.o: modules.o
    265285local_tremain_ini.o: modules.o
     286lpm_advec.o: modules.o
     287lpm_boundary_conds.o: modules.o
     288lpm_calc_liquid_water_content.o: modules.o
    266289lpm_collision_kernels.o: modules.o user_module.o
     290lpm_data_output_particles.o: modules.o
     291lpm_droplet_collision.o: modules.o lpm_collision_kernels.o
     292lpm_droplet_condensation.o: modules.o lpm_collision_kernels.o
     293lpm_exchange_horiz.o: modules.o
     294lpm_extend_particle_array.o: modules.o
     295lpm_extend_tails.o: modules.o
     296lpm_extend_tail_array.o: modules.o
     297lpm_init.o: modules.o lpm_collision_kernels.o random_function.o
     298lpm_init_sgs_tke.o: modules.o
     299lpm_pack_arrays.o: modules.o
     300lpm_read_restart_file.o: modules.o
     301lpm_release_set.o: modules.o random_function.o
     302lpm_set_attributes.o: modules.o
     303lpm_sort_arrays.o: modules.o
     304lpm_write_exchange_statistics.o: modules.o
     305lpm_write_restart_file.o: modules.o
    267306message.o: modules.o
    268307modules.o: modules.f90
     
    271310palm.o: modules.o
    272311parin.o: modules.o
    273 particle_boundary_conds.o: modules.o
    274312plant_canopy_model.o: modules.o
    275313poisfft.o: modules.o fft_xy.o
     
    292330read_var_list.o: modules.o
    293331run_control.o: modules.o
    294 set_particle_attributes.o: modules.o
    295332set_slicer_attributes_dvrp.o: modules.o
    296333singleton.o: singleton.f90
     
    312349user_actions.o: modules.o user_module.o
    313350user_additional_routines.o: modules.o user_module.o
    314 user_advec_particles.o: modules.o user_module.o
    315351user_check_data_output.o: modules.o user_module.o
    316352user_check_data_output_pr.o: modules.o user_module.o
     
    326362user_init_3d_model.o: modules.o user_module.o
    327363user_init_grid.o: modules.o user_module.o
    328 user_init_particles.o: modules.o user_module.o
    329364user_init_plant_canopy.o: modules.o user_module.o
    330365user_last_actions.o: modules.o user_module.o
     366user_lpm_advec.o: modules.o user_module.o
     367user_lpm_init.o: modules.o user_module.o
     368user_lpm_set_attributes.o: modules.o user_module.o
    331369user_module.o: user_module.f90
    332370user_parin.o: modules.o user_module.o
    333 user_particle_attributes.o: modules.o user_module.o
    334371user_read_restart_data.o: modules.o user_module.o
    335372user_spectra.o: modules.o user_module.o
  • palm/trunk/SOURCE/check_open.f90

    r810 r849  
    727727!
    728728!--          Attention: change version number whenever the output format on
    729 !--                     unit 85 is changed (see also in routine advec_particles)
     729!--                     unit 85 is changed (see also in routine
     730!--                     lpm_data_output_particles)
    730731             rtext = 'data format version 3.0'
    731732             WRITE ( 85 )  rtext
  • palm/trunk/SOURCE/init_3d_model.f90

    r826 r849  
    77! Current revisions:
    88! ------------------
    9 !
     9! init_particles renamed lpm_init
    1010!
    1111! Former revisions:
     
    14371437!
    14381438!--    Initialize quantities for handling cloud physics
    1439 !--    This routine must be called before init_particles, because
     1439!--    This routine must be called before lpm_init, because
    14401440!--    otherwise, array pt_d_t, needed in data_output_dvrp (called by
    1441 !--    init_particles) is not defined.
     1441!--    lpm_init) is not defined.
    14421442       CALL init_cloud_physics
    14431443    ENDIF
     
    14451445!
    14461446!-- If required, initialize particles
    1447     IF ( particle_advection )  CALL init_particles
     1447    IF ( particle_advection )  CALL lpm_init
    14481448
    14491449!
  • palm/trunk/SOURCE/lpm.f90

    r848 r849  
    1  SUBROUTINE advec_particles
     1 SUBROUTINE lpm
    22
    33!------------------------------------------------------------------------------!
    44! Current revisions:
    55! ------------------
    6 !
     6! Original routine advec_particles split into several subroutines and renamed
     7! lpm
    78!
    89! Former revisions:
    910! -----------------
    1011! $Id$
    11 !
    12 ! 835 2012-02-22 11:21:19Z raasch
    13 ! Bugfix: array diss can be used only in case of Wang kernel
    1412!
    1513! 831 2012-02-22 00:29:39Z raasch
     
    4341! Splitting of parallel I/O (routine write_particles)
    4442!
    45 ! 667 2010-12-23 12:06:00Z suehring/gryschka
    46 ! Declaration of de_dx, de_dy, de_dz adapted to additional
    47 ! ghost points. Furthermore the calls of exchange_horiz were modified.
    48 !
    49 ! 622 2010-12-10 08:08:13Z raasch
    50 ! optional barriers included in order to speed up collective operations
    51 !
    52 ! 519 2010-03-19 05:30:02Z raasch
    53 ! NetCDF4 output format allows size of particle array to be extended
    54 !
    55 ! 420 2010-01-13 15:10:53Z franke
    56 ! Own weighting factor for every cloud droplet is implemented and condensation
    57 ! and collision of cloud droplets are adjusted accordingly. +delta_v, -s_r3,
    58 ! -s_r4
    59 ! Initialization of variables for the (sub-) timestep is moved to the beginning
    60 ! in order to enable deletion of cloud droplets due to collision processes.
    61 ! Collision efficiency for large cloud droplets has changed according to table
    62 ! of Rogers and Yau. (collision_efficiency)
    63 ! Bugfix: calculation of cloud droplet velocity
    64 ! Bugfix: transfer of particles at south/left edge when new position
    65 !         y>=(ny+0.5)*dy-1.e-12 or x>=(nx+0.5)*dx-1.e-12, very rare
    66 ! Bugfix: calculation of y (collision_efficiency)
    67 !
    68 ! 336 2009-06-10 11:19:35Z raasch
    69 ! Particle attributes are set with new routine set_particle_attributes.
    70 ! Vertical particle advection depends on the particle group.
    71 ! Output of NetCDF messages with aid of message handling routine.
    72 ! Output of messages replaced by message handling routine
    73 ! Bugfix: error in check, if particles moved further than one subdomain length.
    74 !         This check must not be applied for newly released particles
    75 ! Bugfix: several tail counters are initialized, particle_tail_coordinates is
    76 ! only written to file if its third index is > 0
    77 !
    78 ! 212 2008-11-11 09:09:24Z raasch
    79 ! Bugfix in calculating k index in case of oceans runs (sort_particles)
    80 !
    81 ! 150 2008-02-29 08:19:58Z raasch
    82 ! Bottom boundary condition and vertical index calculations adjusted for
    83 ! ocean runs.
    84 !
    85 ! 119 2007-10-17 10:27:13Z raasch
    86 ! Sorting of particles is controlled by dt_sort_particles and moved from
    87 ! the SGS timestep loop after the end of this loop.
    88 ! Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
    89 ! numbers along y
    90 ! Small bugfixes in the SGS part
    91 !
    92 ! 106 2007-08-16 14:30:26Z raasch
    93 ! remaining variables iran changed to iran_part
    94 !
    95 ! 95 2007-06-02 16:48:38Z raasch
    96 ! hydro_press renamed hyp
    97 !
    98 ! 75 2007-03-22 09:54:05Z raasch
    99 ! Particle reflection at vertical walls implemented in new subroutine
    100 ! particle_boundary_conds,
    101 ! vertical walls are regarded in the SGS model,
    102 ! + user_advec_particles, particles-package is now part of the defaut code,
    103 ! array arguments in sendrecv calls have to refer to first element (1) due to
    104 ! mpich (mpiI) interface requirements,
    105 ! 2nd+3rd argument removed from exchange horiz
    106 !
    107 ! 16 2007-02-15 13:16:47Z raasch
    108 ! Bugfix: wrong if-clause from revision 1.32
    109 !
    110 ! r4 | raasch | 2007-02-13 12:33:16 +0100 (Tue, 13 Feb 2007)
    111 ! RCS Log replace by Id keyword, revision history cleaned up
    112 !
    113 ! Revision 1.32  2007/02/11 12:48:20  raasch
    114 ! Allways the lower level k is used for interpolation
    115 ! Bugfix: new particles are released only if end_time_prel > simulated_time
    116 ! Bugfix: transfer of particles when x < -0.5*dx (0.0 before), etc.,
    117 !         index i,j used instead of cartesian (x,y) coordinate to check for
    118 !         transfer because this failed under very rare conditions
    119 ! Bugfix: calculation of number of particles with same radius as the current
    120 !         particle (cloud droplet code)
    121 !
    122 ! Revision 1.31  2006/08/17 09:21:01  raasch
    123 ! Two more compilation errors removed from the last revision
    124 !
    125 ! Revision 1.30  2006/08/17 09:11:17  raasch
    126 ! Two compilation errors removed from the last revision
    127 !
    128 ! Revision 1.29  2006/08/04 14:05:01  raasch
    129 ! Subgrid scale velocities are (optionally) included for calculating the
    130 ! particle advection, new counters trlp_count_sum, etc. for accumulating
    131 ! the number of particles exchanged between the subdomains during all
    132 ! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z,
    133 ! izuf renamed iran, output of particle time series
    134 !
    13543! Revision 1.1  1999/11/25 16:16:06  raasch
    13644! Initial revision
     
    14351
    14452    USE arrays_3d
    145     USE cloud_parameters
    146     USE constants
    14753    USE control_parameters
    14854    USE cpulog
    149     USE grid_variables
    150     USE indices
    15155    USE interfaces
    152     USE lpm_collision_kernels_mod
    153     USE netcdf_control
    15456    USE particle_attributes
    15557    USE pegrid
    156     USE random_function_mod
    15758    USE statistics
    15859
    15960    IMPLICIT NONE
    16061
    161     INTEGER ::  agp, deleted_particles, deleted_tails, eclass, i, ie, ii, inc, &
    162                 internal_timestep_count, is, j, jj, js, jtry, k, kk, kw, m, n, &
    163                 nc, nd, nn, num_gp, pse, psi, rclass_l, rclass_s,  &
    164                 tlength, trlp_count, trlp_count_sum, trlp_count_recv,          &
    165                 trlp_count_recv_sum, trlpt_count, trlpt_count_recv,            &
    166                 trnp_count, trnp_count_sum, trnp_count_recv,                   &
    167                 trnp_count_recv_sum, trnpt_count, trnpt_count_recv,            &
    168                 trrp_count, trrp_count_sum, trrp_count_recv,                   &
    169                 trrp_count_recv_sum, trrpt_count, trrpt_count_recv,            &
    170                 trsp_count, trsp_count_sum, trsp_count_recv,                   &
    171                 trsp_count_recv_sum, trspt_count, trspt_count_recv
    172 
    173     INTEGER ::  gp_outside_of_building(1:8)
    174 
    175     INTEGER, PARAMETER ::  maxtry = 40   ! for Rosenbrock method
    176 
    177     LOGICAL ::  dt_3d_reached, dt_3d_reached_l, prt_position
    178 
    179     REAL    ::  aa, afactor, arg, bb, cc, dd, ddenom, delta_r, delta_v, &
    180                 dens_ratio, de_dt, de_dt_min, de_dx_int, de_dx_int_l, &
    181                 de_dx_int_u, de_dy_int, de_dy_int_l, de_dy_int_u, de_dz_int, &
    182                 de_dz_int_l, de_dz_int_u, diss_int, diss_int_l, diss_int_u, &
    183                 distance, drdt, drdt_ini, drdt_m, dt_ros, dt_ros_last, &
    184                 dt_ros_next, dt_ros_sum, dt_ros_sum_ini, d2rdt2, d2rdtdr, &
    185                 dt_gap, dt_particle, dt_particle_m, d_sum, epsilon, e_a, e_int,&
    186                 e_int_l, e_int_u, e_mean_int, e_s, err_ros, errmax, exp_arg, &
    187                 exp_term, fs_int, gg, g1, g2, g3, g4, integral, &
    188                 lagr_timescale, lw_max, mean_r, new_r, p_int, pt_int, &
    189                 pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, &
    190                 ql_int_u, random_gauss, r_ros, r_ros_ini, &
    191                 sigma, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,vv_int, &
    192                 v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u, x, y
    193 !
    194 !-- Parameters for Rosenbrock method
    195     REAL, PARAMETER ::  a21 = 2.0, a31 = 48.0/25.0, a32 = 6.0/25.0,        &
    196                         a2x = 1.0, a3x = 3.0/5.0, b1 = 19.0/9.0, b2 = 0.5, &
    197                         b3 = 25.0/108.0, b4 = 125.0/108.0, c21 = -8.0,     &
    198                         c31 = 372.0/25.0, c32 = 12.0/5.0,                  &
    199                         c41 = -112.0/125.0, c42 = -54.0/125.0,             &
    200                         c43 = -2.0/5.0, c1x = 0.5, c2x= -3.0/2.0,          &
    201                         c3x = 121.0/50.0, c4x = 29.0/250.0,                &
    202                         errcon = 0.1296, e1 = 17.0/54.0, e2 = 7.0/36.0,    &
    203                         e3 = 0.0, e4 = 125.0/108.0, gam = 0.5, grow = 1.5, &
    204                         pgrow = -0.25, pshrnk = -1.0/3.0, shrnk = 0.5
    205 
    206     REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
    207 
    208     REAL    ::  location(1:30,1:3)
    209 
    210     REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  de_dx, de_dy, de_dz
    211 
    212     REAL, DIMENSION(:,:,:), ALLOCATABLE ::  trlpt, trnpt, trrpt, trspt
    213 
    214     TYPE(particle_type) ::  tmp_particle
    215 
    216     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp, trnp, trrp, trsp
    217 
    218 
    219     CALL cpu_log( log_point(25), 'advec_particles', 'start' )
    220 
    221 
    222 !
    223 !-- Write particle data on file for later analysis.
    224 !-- This has to be done here (before particles are advected) in order
    225 !-- to allow correct output in case of dt_write_particle_data = dt_prel =
    226 !-- particle_maximum_age. Otherwise (if output is done at the end of this
    227 !-- subroutine), the relevant particles would have been already deleted.
    228 !-- The MOD function allows for changes in the output interval with restart
    229 !-- runs.
    230 !-- Attention: change version number for unit 85 (in routine check_open)
    231 !--            whenever the output format for this unit is changed!
     62    INTEGER ::  m
     63
     64
     65    CALL cpu_log( log_point(25), 'lpm', 'start' )
     66
     67!
     68!-- Write particle data at current time on file.
     69!-- This has to be done here, before particles are further processed,
     70!-- because they may be deleted within this timestep (in case that
     71!-- dt_write_particle_data = dt_prel = particle_maximum_age).
    23272    time_write_particle_data = time_write_particle_data + dt_3d
    23373    IF ( time_write_particle_data >= dt_write_particle_data )  THEN
    23474
    235        CALL cpu_log( log_point_s(40), 'advec_part_io', 'start' )
    236        CALL check_open( 85 )
    237        WRITE ( 85 )  simulated_time, maximum_number_of_particles, &
    238                      number_of_particles
    239        WRITE ( 85 )  particles
    240        WRITE ( 85 )  maximum_number_of_tailpoints, maximum_number_of_tails, &
    241                      number_of_tails
    242        IF ( maximum_number_of_tails > 0 )  THEN
    243           WRITE ( 85 )  particle_tail_coordinates
    244        ENDIF
    245        CALL close_file( 85 )
    246 
    247        IF ( netcdf_output )  CALL output_particles_netcdf
    248 
     75       CALL lpm_data_output_particles
     76!
     77!--    The MOD function allows for changes in the output interval with restart
     78!--    runs.
    24979       time_write_particle_data = MOD( time_write_particle_data, &
    25080                                  MAX( dt_write_particle_data, dt_3d ) )
    251        CALL cpu_log( log_point_s(40), 'advec_part_io', 'stop' )
    252     ENDIF
    253 
    254 !
    255 !--    Initialize variables for the (sub-) timestep, i.e. for
    256 !--    marking those particles to be deleted after the timestep
    257        particle_mask     = .TRUE.
    258        deleted_particles = 0
    259        trlp_count_recv   = 0
    260        trnp_count_recv   = 0
    261        trrp_count_recv   = 0
    262        trsp_count_recv   = 0
    263        trlpt_count_recv  = 0
    264        trnpt_count_recv  = 0
    265        trrpt_count_recv  = 0
    266        trspt_count_recv  = 0
    267        IF ( use_particle_tails )  THEN
    268           tail_mask     = .TRUE.
    269        ENDIF
    270        deleted_tails = 0
    271 
    272 !
    273 !-- Calculate exponential term used in case of particle inertia for each
    274 !-- of the particle groups
    275     CALL cpu_log( log_point_s(41), 'advec_part_exp', 'start' )
    276     DO  m = 1, number_of_particle_groups
    277        IF ( particle_groups(m)%density_ratio /= 0.0 )  THEN
    278           particle_groups(m)%exp_arg  =                                        &
    279                     4.5 * particle_groups(m)%density_ratio *                   &
    280                     molecular_viscosity / ( particle_groups(m)%radius )**2
    281           particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * &
    282                                              dt_3d )
    283        ENDIF
    284     ENDDO
    285     CALL cpu_log( log_point_s(41), 'advec_part_exp', 'stop' )
    286 
    287 !
    288 !-- Particle (droplet) growth by condensation/evaporation and collision
    289     IF ( cloud_droplets )  THEN
    290 
    291 !
    292 !--    Reset summation arrays
    293        ql_c = 0.0;  ql_v = 0.0;  ql_vp = 0.0
    294        delta_r=0.0;  delta_v=0.0
    295 
    296 !
    297 !--    Particle growth by condensation/evaporation
    298        CALL cpu_log( log_point_s(42), 'advec_part_cond', 'start' )
    299        DO  n = 1, number_of_particles
    300 !
    301 !--       Interpolate temperature and humidity.
    302 !--       First determine left, south, and bottom index of the arrays.
    303           i = particles(n)%x * ddx
    304           j = particles(n)%y * ddy
    305           k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    306               + offset_ocean_nzt                     ! only exact if equidistant
    307 
    308           x  = particles(n)%x - i * dx
    309           y  = particles(n)%y - j * dy
    310           aa = x**2          + y**2
    311           bb = ( dx - x )**2 + y**2
    312           cc = x**2          + ( dy - y )**2
    313           dd = ( dx - x )**2 + ( dy - y )**2
    314           gg = aa + bb + cc + dd
    315 
    316           pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
    317                      + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
    318                      ) / ( 3.0 * gg )
    319 
    320           pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
    321                      + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
    322                      ) / ( 3.0 * gg )
    323 
    324           pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
    325                               ( pt_int_u - pt_int_l )
    326 
    327           q_int_l = ( ( gg - aa ) * q(k,j,i)   + ( gg - bb ) * q(k,j,i+1)   &
    328                     + ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) &
    329                     ) / ( 3.0 * gg )
    330 
    331           q_int_u = ( ( gg-aa ) * q(k+1,j,i)   + ( gg-bb ) * q(k+1,j,i+1)   &
    332                     + ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) &
    333                     ) / ( 3.0 * gg )
    334 
    335           q_int = q_int_l + ( particles(n)%z - zu(k) ) / dz * &
    336                               ( q_int_u - q_int_l )
    337 
    338           ql_int_l = ( ( gg - aa ) * ql(k,j,i)   + ( gg - bb ) * ql(k,j,i+1)   &
    339                      + ( gg - cc ) * ql(k,j+1,i) + ( gg - dd ) * ql(k,j+1,i+1) &
    340                      ) / ( 3.0 * gg )
    341 
    342           ql_int_u = ( ( gg-aa ) * ql(k+1,j,i)   + ( gg-bb ) * ql(k+1,j,i+1)   &
    343                      + ( gg-cc ) * ql(k+1,j+1,i) + ( gg-dd ) * ql(k+1,j+1,i+1) &
    344                      ) / ( 3.0 * gg )
    345 
    346           ql_int = ql_int_l + ( particles(n)%z - zu(k) ) / dz * &
    347                                 ( ql_int_u - ql_int_l )
    348 
    349 !
    350 !--       Calculate real temperature and saturation vapor pressure
    351           p_int = hyp(k) + ( particles(n)%z - zu(k) ) / dz * ( hyp(k+1)-hyp(k) )
    352           t_int = pt_int * ( p_int / 100000.0 )**0.286
    353 
    354           e_s = 611.0 * EXP( l_d_rv * ( 3.6609E-3 - 1.0 / t_int ) )
    355 
    356 !
    357 !--       Current vapor pressure
    358           e_a = q_int * p_int / ( 0.378 * q_int + 0.622 )
    359 
    360 !
    361 !--       Thermal conductivity for water (from Rogers and Yau, Table 7.1),
    362 !--       diffusivity for water vapor (after Hall und Pruppacher, 1976)
    363           thermal_conductivity_l = 7.94048E-05 * t_int + 0.00227011
    364           diff_coeff_l           = 0.211E-4 * ( t_int / 273.15 )**1.94 * &
    365                                    ( 101325.0 / p_int)
    366 !
    367 !--       Change in radius by condensation/evaporation
    368           IF ( particles(n)%radius >= 1.0E-6  .OR.  &
    369                .NOT. curvature_solution_effects )  THEN
    370 !
    371 !--          Approximation for large radii, where curvature and solution
    372 !--          effects can be neglected
    373              arg = particles(n)%radius**2 + 2.0 * dt_3d *                     &
    374                            ( e_a / e_s - 1.0 ) /                              &
    375                            ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
    376                              thermal_conductivity_l +                         &
    377                              rho_l * r_v * t_int / diff_coeff_l / e_s )
    378              IF ( arg < 1.0E-16 )  THEN
    379                 new_r = 1.0E-8
    380              ELSE
    381                 new_r = SQRT( arg )
    382              ENDIF
    383           ENDIF
    384 
    385           IF ( curvature_solution_effects  .AND.                              &
    386                ( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) &
    387           THEN
    388 !
    389 !--          Curvature and solutions effects are included in growth equation.
    390 !--          Change in Radius is calculated with 4th-order Rosenbrock method
    391 !--          for stiff o.d.e's with monitoring local truncation error to adjust
    392 !--          stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
    393 !--          For larger radii the simple analytic method (see ELSE) gives
    394 !--          almost the same results.
    395 !
    396 !--          Surface tension after (Straka, 2009)
    397              sigma = 0.0761 - 0.00155 * ( t_int - 273.15 )
    398 
    399              r_ros = particles(n)%radius
    400              dt_ros_sum  = 0.0      ! internal integrated time (s)
    401              internal_timestep_count = 0
    402 
    403              ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
    404                                ( l_v / ( r_v * t_int ) - 1.0 ) *               &
    405                                rho_l * l_v / ( thermal_conductivity_l * t_int )&
    406                              )
    407              
    408              afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
    409 
    410              IF ( particles(n)%rvar3 == -9999999.9 )  THEN
    411 !
    412 !--             First particle timestep. Derivative has to be calculated.
    413                 drdt_m  = ddenom / r_ros * ( e_a / e_s - 1.0 -    &
    414                                              afactor / r_ros +    &
    415                                              bfactor / r_ros**3 )
    416              ELSE
    417 !
    418 !--             Take value from last PALM timestep
    419                 drdt_m = particles(n)%rvar3
    420              ENDIF
    421 !
    422 !--          Take internal timestep values from the end of last PALM timestep
    423              dt_ros_last = particles(n)%rvar1
    424              dt_ros_next = particles(n)%rvar2
    425 !
    426 !--          Internal timestep must not be larger than PALM timestep
    427              dt_ros = MIN( dt_ros_next, dt_3d )
    428 !
    429 !--          Integrate growth equation in time unless PALM timestep is reached
    430              DO WHILE ( dt_ros_sum < dt_3d )
    431 
    432                 internal_timestep_count = internal_timestep_count + 1
    433 
    434 !
    435 !--             Derivative at starting value
    436                 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
    437                                           bfactor / r_ros**3 )
    438                 drdt_ini       = drdt
    439                 dt_ros_sum_ini = dt_ros_sum
    440                 r_ros_ini      = r_ros
    441 
    442 !
    443 !--             Calculate time derivative of dr/dt
    444                 d2rdt2 = ( drdt - drdt_m ) / dt_ros_last
    445 
    446 !
    447 !--             Calculate radial derivative of dr/dt
    448                 d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
    449                                      2.0 * afactor / r_ros**3 -     &
    450                                      4.0 * bfactor / r_ros**5 )
    451 !
    452 !--             Adjust stepsize unless required accuracy is reached
    453                 DO  jtry = 1, maxtry+1
    454 
    455                    IF ( jtry == maxtry+1 )  THEN
    456                       message_string = 'maxtry > 40 in Rosenbrock method'
    457                       CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 )
    458                    ENDIF
    459 
    460                    aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
    461                    g1    = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa
    462                    r_ros = r_ros_ini + a21 * g1
    463                    drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
    464                                               afactor / r_ros + &
    465                                               bfactor / r_ros**3 )
    466 
    467                    g2    = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )&
    468                            / aa
    469                    r_ros = r_ros_ini + a31 * g1 + a32 * g2
    470                    drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
    471                                               afactor / r_ros + &
    472                                               bfactor / r_ros**3 )
    473 
    474                    g3    = ( drdt + dt_ros * c3x * d2rdt2 + &
    475                              ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
    476                    g4    = ( drdt + dt_ros * c4x * d2rdt2 + &
    477                              ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
    478                    r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
    479 
    480                    dt_ros_sum = dt_ros_sum_ini + dt_ros
    481 
    482                    IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
    483                       message_string = 'zero stepsize in Rosenbrock method'
    484                       CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 )
    485                    ENDIF
    486 !
    487 !--                Calculate error
    488                    err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4
    489                    errmax  = 0.0
    490                    errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
    491 !
    492 !--                Leave loop if accuracy is sufficient, otherwise try again
    493 !--                with a reduced stepsize
    494                    IF ( errmax <= 1.0 )  THEN
    495                       EXIT
    496                    ELSE
    497                       dt_ros = SIGN(                                           &
    498                                     MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
    499                                          shrnk * ABS( dt_ros ) ),              &
    500                                     dt_ros                                     &
    501                                    )
    502                    ENDIF
    503 
    504                 ENDDO  ! loop for stepsize adjustment
    505 
    506 !
    507 !--             Calculate next internal timestep
    508                 IF ( errmax > errcon )  THEN
    509                    dt_ros_next = 0.9 * dt_ros * errmax**pgrow
    510                 ELSE
    511                    dt_ros_next = grow * dt_ros
    512                 ENDIF
    513 
    514 !
    515 !--             Estimated timestep is reduced if the PALM time step is exceeded
    516                 dt_ros_last = dt_ros
    517                 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
    518                    dt_ros = dt_3d - dt_ros_sum
    519                 ELSE
    520                    dt_ros = dt_ros_next
    521                 ENDIF
    522 
    523                 drdt_m = drdt
    524 
    525              ENDDO
    526 !
    527 !--          Store derivative and internal timestep values for next PALM step
    528              particles(n)%rvar1 = dt_ros_last
    529              particles(n)%rvar2 = dt_ros_next
    530              particles(n)%rvar3 = drdt_m
    531 
    532              new_r = r_ros
    533 !
    534 !--          Radius should not fall below 1E-8 because Rosenbrock method may
    535 !--          lead to errors otherwise
    536              new_r = MAX( new_r, 1.0E-8 )
    537 
    538           ENDIF
    539 
    540           delta_r = new_r - particles(n)%radius
    541 
    542 !
    543 !--       Sum up the change in volume of liquid water for the respective grid
    544 !--       volume (this is needed later on for calculating the release of
    545 !--       latent heat)
    546           i = ( particles(n)%x + 0.5 * dx ) * ddx
    547           j = ( particles(n)%y + 0.5 * dy ) * ddy
    548           k = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
    549               ! only exact if equidistant
    550 
    551           ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor *            &
    552                                       rho_l * 1.33333333 * pi *               &
    553                                       ( new_r**3 - particles(n)%radius**3 ) / &
    554                                       ( rho_surface * dx * dy * dz )
    555           IF ( ql_c(k,j,i) > 100.0 )  THEN
    556              WRITE( message_string, * ) 'k=',k,' j=',j,' i=',i,      &
    557                           ' ql_c=',ql_c(k,j,i), ' &part(',n,')%wf=', &
    558                           particles(n)%weight_factor,' delta_r=',delta_r
    559              CALL message( 'advec_particles', 'PA0143', 2, 2, -1, 6, 1 )
    560           ENDIF
    561 
    562 !
    563 !--       Change the droplet radius
    564           IF ( ( new_r - particles(n)%radius ) < 0.0  .AND.  new_r < 0.0 ) &
    565           THEN
    566              WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i,   &
    567                           ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,  &
    568                           ' &delta_r=',delta_r,                      &
    569                           ' particle_radius=',particles(n)%radius
    570              CALL message( 'advec_particles', 'PA0144', 2, 2, -1, 6, 1 )
    571           ENDIF
    572 
    573 !
    574 !--       Sum up the total volume of liquid water (needed below for
    575 !--       re-calculating the weighting factors)
    576           ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r**3
    577 
    578           particles(n)%radius = new_r
    579 
    580 !
    581 !--       Determine radius class of the particle needed for collision
    582           IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  use_kernel_tables ) &
    583           THEN
    584              particles(n)%class = ( LOG( new_r ) - rclass_lbound ) /  &
    585                                   ( rclass_ubound - rclass_lbound ) * &
    586                                   radius_classes
    587              particles(n)%class = MIN( particles(n)%class, radius_classes )
    588              particles(n)%class = MAX( particles(n)%class, 1 )
    589           ENDIF
    590 
    591        ENDDO
    592        CALL cpu_log( log_point_s(42), 'advec_part_cond', 'stop' )
    593 
    594 !
    595 !--    Particle growth by collision
    596        IF ( collision_kernel /= 'none' )  THEN
    597 
    598        CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
    599 
    600        DO  i = nxl, nxr
    601           DO  j = nys, nyn
    602              DO  k = nzb+1, nzt
    603 !
    604 !--             Collision requires at least two particles in the box
    605                 IF ( prt_count(k,j,i) > 1 )  THEN
    606 !
    607 !--                First, sort particles within the gridbox by their size,
    608 !--                using Shell's method (see Numerical Recipes)
    609 !--                NOTE: In case of using particle tails, the re-sorting of
    610 !--                ----  tails would have to be included here!
    611                    psi = prt_start_index(k,j,i) - 1
    612                    inc = 1
    613                    DO WHILE ( inc <= prt_count(k,j,i) )
    614                       inc = 3 * inc + 1
    615                    ENDDO
    616 
    617                    DO WHILE ( inc > 1 )
    618                       inc = inc / 3
    619                       DO  is = inc+1, prt_count(k,j,i)
    620                          tmp_particle = particles(psi+is)
    621                          js = is
    622                          DO WHILE ( particles(psi+js-inc)%radius >             &
    623                                     tmp_particle%radius )
    624                             particles(psi+js) = particles(psi+js-inc)
    625                             js = js - inc
    626                             IF ( js <= inc )  EXIT
    627                          ENDDO
    628                          particles(psi+js) = tmp_particle
    629                       ENDDO
    630                    ENDDO
    631 
    632                    psi = prt_start_index(k,j,i)
    633                    pse = psi + prt_count(k,j,i)-1
    634 
    635 !
    636 !--                Now apply the different kernels
    637                    IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
    638                         use_kernel_tables )  THEN
    639 !
    640 !--                   Fast method with pre-calculated efficiencies for
    641 !--                   discrete radius- and dissipation-classes.
    642 !
    643 !--                   Determine dissipation class index of this gridbox
    644                       IF ( wang_kernel )  THEN
    645                          eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * &
    646                                        dissipation_classes ) + 1
    647                          epsilon = diss(k,j,i)
    648                       ELSE
    649                          epsilon = 0.0
    650                       ENDIF
    651                       IF ( hall_kernel  .OR.  epsilon * 1.0E4 < 0.001 )  THEN
    652                          eclass = 0   ! Hall kernel is used
    653                       ELSE
    654                          eclass = MIN( dissipation_classes, eclass )
    655                       ENDIF
    656 
    657                       DO  n = pse, psi+1, -1
    658 
    659                          integral = 0.0
    660                          lw_max   = 0.0
    661                          rclass_l = particles(n)%class
    662 !
    663 !--                      Calculate growth of collector particle radius using all
    664 !--                      droplets smaller than current droplet
    665                          DO  is = psi, n-1
    666 
    667                             rclass_s = particles(is)%class
    668                             integral = integral +                             &
    669                                        ( particles(is)%radius**3 *            &
    670                                          ckernel(rclass_l,rclass_s,eclass) *  &
    671                                          particles(is)%weight_factor )
    672 !
    673 !--                         Calculate volume of liquid water of the collected
    674 !--                         droplets which is the maximum liquid water available
    675 !--                         for droplet growth   
    676                             lw_max = lw_max + ( particles(is)%radius**3 * &
    677                                                 particles(is)%weight_factor )
    678 
    679                          ENDDO
    680 
    681 !
    682 !--                      Change in radius of collector droplet due to collision
    683                          delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
    684                                    integral * dt_3d * ddx * ddy / dz
    685 
    686 !
    687 !--                      Change in volume of collector droplet due to collision
    688                          delta_v = particles(n)%weight_factor                  &
    689                                    * ( ( particles(n)%radius + delta_r )**3    &
    690                                        - particles(n)%radius**3 )
    691 
    692                          IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
    693 !-- replace by message call
    694                             PRINT*, 'Particle has grown to much because the',  &
    695                                     ' change of volume of particle is larger', &
    696                                     ' than liquid water available!'
    697 
    698                          ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
    699 !-- can this case really happen??
    700                             DO  is = psi, n-1
    701 
    702                                particles(is)%weight_factor = 0.0
    703                                particle_mask(is)  = .FALSE.
    704                                deleted_particles = deleted_particles + 1
    705 
    706                             ENDDO
    707 
    708                          ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
    709 !
    710 !--                         Calculate new weighting factor of collected droplets
    711                             DO  is = psi, n-1
    712 
    713                                rclass_s = particles(is)%class
    714                                particles(is)%weight_factor = &
    715                                   particles(is)%weight_factor - &
    716                                   ( ( ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor &
    717                                     / integral ) * delta_v )
    718 
    719                                IF ( particles(is)%weight_factor < 0.0 )  THEN
    720                                      WRITE( message_string, * ) 'negative ',   &
    721                                                         'weighting factor: ',  &
    722                                                      particles(is)%weight_factor
    723                                      CALL message( 'advec_particles', '', 2, 2,&
    724                                                    -1, 6, 1 )
    725 
    726                                ELSEIF ( particles(is)%weight_factor == 0.0 )  &
    727                                THEN
    728 
    729                                    particles(is)%weight_factor = 0.0
    730                                    particle_mask(is)  = .FALSE.
    731                                    deleted_particles = deleted_particles + 1
    732 
    733                                ENDIF
    734 
    735                             ENDDO
    736 
    737                          ENDIF
    738 
    739                          particles(n)%radius = particles(n)%radius + delta_r
    740                          ql_vp(k,j,i) = ql_vp(k,j,i) + &
    741                                                     particles(n)%weight_factor &
    742                                                     * particles(n)%radius**3
    743 
    744                       ENDDO
    745 
    746                    ELSEIF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
    747                             .NOT. use_kernel_tables )  THEN
    748 !
    749 !--                   Collision efficiencies are calculated for every new
    750 !--                   grid box. First, allocate memory for kernel table.
    751 !--                   Third dimension is 1, because table is re-calculated for
    752 !--                   every new dissipation value.
    753                       ALLOCATE( ckernel(prt_start_index(k,j,i):                &
    754                                prt_start_index(k,j,i)+prt_count(k,j,i)-1,      &
    755                                prt_start_index(k,j,i):                         &
    756                                prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) )
    757 !
    758 !--                   Now calculate collision efficiencies for this box
    759                       CALL recalculate_kernel( i, j, k )
    760 
    761                       DO  n = pse, psi+1, -1
    762 
    763                          integral = 0.0
    764                          lw_max   = 0.0
    765 !
    766 !--                      Calculate growth of collector particle radius using all
    767 !--                      droplets smaller than current droplet
    768                          DO  is = psi, n-1
    769 
    770                             integral = integral + particles(is)%radius**3 * &
    771                                                   ckernel(n,is,1) *         &
    772                                                   particles(is)%weight_factor
    773 !
    774 !--                         Calculate volume of liquid water of the collected
    775 !--                         droplets which is the maximum liquid water available
    776 !--                         for droplet growth   
    777                             lw_max = lw_max + ( particles(is)%radius**3 * &
    778                                                 particles(is)%weight_factor )
    779 
    780                          ENDDO
    781 
    782 !
    783 !--                      Change in radius of collector droplet due to collision
    784                          delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
    785                                    integral * dt_3d * ddx * ddy / dz
    786 
    787 !
    788 !--                      Change in volume of collector droplet due to collision
    789                          delta_v = particles(n)%weight_factor                  &
    790                                    * ( ( particles(n)%radius + delta_r )**3    &
    791                                        - particles(n)%radius**3 )
    792 
    793                          IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
    794 !-- replace by message call
    795                             PRINT*, 'Particle has grown to much because the',  &
    796                                     ' change of volume of particle is larger', &
    797                                     ' than liquid water available!'
    798 
    799                          ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
    800 !-- can this case really happen??
    801                             DO  is = psi, n-1
    802 
    803                                particles(is)%weight_factor = 0.0
    804                                particle_mask(is)  = .FALSE.
    805                                deleted_particles = deleted_particles + 1
    806 
    807                             ENDDO
    808 
    809                          ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
    810 !
    811 !--                         Calculate new weighting factor of collected droplets
    812                             DO  is = psi, n-1
    813 
    814                                particles(is)%weight_factor =                   &
    815                                       particles(is)%weight_factor -            &
    816                                       ( ckernel(n,is,1) / integral * delta_v * &
    817                                         particles(is)%weight_factor )
    818 
    819                                IF ( particles(is)%weight_factor < 0.0 )  THEN
    820                                      WRITE( message_string, * ) 'negative ',   &
    821                                                         'weighting factor: ',  &
    822                                                      particles(is)%weight_factor
    823                                      CALL message( 'advec_particles', '', 2, 2,&
    824                                                    -1, 6, 1 )
    825 
    826                                ELSEIF ( particles(is)%weight_factor == 0.0 )  &
    827                                THEN
    828 
    829                                    particles(is)%weight_factor = 0.0
    830                                    particle_mask(is)  = .FALSE.
    831                                    deleted_particles = deleted_particles + 1
    832 
    833                                ENDIF
    834 
    835                             ENDDO
    836 
    837                          ENDIF
    838 
    839                          particles(n)%radius = particles(n)%radius + delta_r
    840                          ql_vp(k,j,i) = ql_vp(k,j,i) + &
    841                                                     particles(n)%weight_factor &
    842                                                     * particles(n)%radius**3
    843 
    844                       ENDDO
    845 
    846                       DEALLOCATE( ckernel )
    847 
    848                    ELSEIF ( palm_kernel )  THEN
    849 !
    850 !--                   PALM collision kernel
    851 !
    852 !--                   Calculate the mean radius of all those particles which
    853 !--                   are of smaller size than the current particle and
    854 !--                   use this radius for calculating the collision efficiency
    855                       DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
    856 
    857                          sl_r3 = 0.0
    858                          sl_r4 = 0.0
    859 
    860                          DO is = n-1, psi, -1
    861                             IF ( particles(is)%radius < particles(n)%radius )  &
    862                             THEN
    863                                sl_r3 = sl_r3 + particles(is)%weight_factor     &
    864                                              *( particles(is)%radius**3 )
    865                                sl_r4 = sl_r4 + particles(is)%weight_factor     &
    866                                              *( particles(is)%radius**4 )
    867                             ENDIF
    868                          ENDDO
    869 
    870                          IF ( ( sl_r3 ) > 0.0 )  THEN
    871                             mean_r = ( sl_r4 ) / ( sl_r3 )
    872 
    873                             CALL collision_efficiency( mean_r,                 &
    874                                                     particles(n)%radius,       &
    875                                                     effective_coll_efficiency )
    876 
    877                          ELSE
    878                             effective_coll_efficiency = 0.0
    879                          ENDIF
    880 
    881                          IF ( effective_coll_efficiency > 1.0  .OR.            &
    882                               effective_coll_efficiency < 0.0 )                &
    883                          THEN
    884                             WRITE( message_string, * )  'collision_efficien' , &
    885                                    'cy out of range:' ,effective_coll_efficiency
    886                             CALL message( 'advec_particles', 'PA0145', 2, 2,   &
    887                                           -1, 6, 1 )
    888                          ENDIF
    889 
    890 !
    891 !--                      Interpolation of ...
    892                          ii = particles(n)%x * ddx
    893                          jj = particles(n)%y * ddy
    894                          kk = ( particles(n)%z + 0.5 * dz ) / dz
    895 
    896                          x  = particles(n)%x - ii * dx
    897                          y  = particles(n)%y - jj * dy
    898                          aa = x**2          + y**2
    899                          bb = ( dx - x )**2 + y**2
    900                          cc = x**2          + ( dy - y )**2
    901                          dd = ( dx - x )**2 + ( dy - y )**2
    902                          gg = aa + bb + cc + dd
    903 
    904                          ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
    905                                                               ql(kk,jj,ii+1)   &
    906                                     + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
    907                                                               ql(kk,jj+1,ii+1) &
    908                                     ) / ( 3.0 * gg )
    909 
    910                          ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
    911                                                             ql(kk+1,jj,ii+1)   &
    912                                     + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
    913                                                             ql(kk+1,jj+1,ii+1) &
    914                                     ) / ( 3.0 * gg )
    915 
    916                          ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
    917                                              ( ql_int_u - ql_int_l )
    918 
    919 !
    920 !--                      Interpolate u velocity-component
    921                          ii = ( particles(n)%x + 0.5 * dx ) * ddx
    922                          jj =   particles(n)%y * ddy
    923                          kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist
    924 
    925                          IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1
    926 
    927                          x  = particles(n)%x + ( 0.5 - ii ) * dx
    928                          y  = particles(n)%y - jj * dy
    929                          aa = x**2          + y**2
    930                          bb = ( dx - x )**2 + y**2
    931                          cc = x**2          + ( dy - y )**2
    932                          dd = ( dx - x )**2 + ( dy - y )**2
    933                          gg = aa + bb + cc + dd
    934 
    935                          u_int_l = ( (gg-aa) * u(kk,jj,ii)   + (gg-bb) *       &
    936                                                               u(kk,jj,ii+1)    &
    937                                    + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
    938                                                               u(kk,jj+1,ii+1)  &
    939                                    ) / ( 3.0 * gg ) - u_gtrans
    940                          IF ( kk+1 == nzt+1 )  THEN
    941                             u_int = u_int_l
    942                          ELSE
    943                             u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
    944                                                              u(kk+1,jj,ii+1)   &
    945                                       + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
    946                                                              u(kk+1,jj+1,ii+1) &
    947                                       ) / ( 3.0 * gg ) - u_gtrans
    948                             u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
    949                                               * ( u_int_u - u_int_l )
    950                          ENDIF
    951 
    952 !
    953 !--                      Same procedure for interpolation of the v velocity-com-
    954 !--                      ponent (adopt index k from u velocity-component)
    955                          ii =   particles(n)%x * ddx
    956                          jj = ( particles(n)%y + 0.5 * dy ) * ddy
    957 
    958                          x  = particles(n)%x - ii * dx
    959                          y  = particles(n)%y + ( 0.5 - jj ) * dy
    960                          aa = x**2          + y**2
    961                          bb = ( dx - x )**2 + y**2
    962                          cc = x**2          + ( dy - y )**2
    963                          dd = ( dx - x )**2 + ( dy - y )**2
    964                          gg = aa + bb + cc + dd
    965 
    966                          v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *   &
    967                                                                v(kk,jj,ii+1)   &
    968                                    + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
    969                                                                v(kk,jj+1,ii+1) &
    970                                    ) / ( 3.0 * gg ) - v_gtrans
    971                          IF ( kk+1 == nzt+1 )  THEN
    972                             v_int = v_int_l
    973                          ELSE
    974                             v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
    975                                                                v(kk+1,jj,ii+1) &
    976                                       + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
    977                                                              v(kk+1,jj+1,ii+1) &
    978                                       ) / ( 3.0 * gg ) - v_gtrans
    979                             v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
    980                                               * ( v_int_u - v_int_l )
    981                          ENDIF
    982 
    983 !
    984 !--                      Same procedure for interpolation of the w velocity-com-
    985 !--                      ponent (adopt index i from v velocity-component)
    986                          jj = particles(n)%y * ddy
    987                          kk = particles(n)%z / dz
    988 
    989                          x  = particles(n)%x - ii * dx
    990                          y  = particles(n)%y - jj * dy
    991                          aa = x**2          + y**2
    992                          bb = ( dx - x )**2 + y**2
    993                          cc = x**2          + ( dy - y )**2
    994                          dd = ( dx - x )**2 + ( dy - y )**2
    995                          gg = aa + bb + cc + dd
    996 
    997                          w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *   &
    998                                                                  w(kk,jj,ii+1) &
    999                                    + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
    1000                                                                w(kk,jj+1,ii+1) &
    1001                                    ) / ( 3.0 * gg )
    1002                          IF ( kk+1 == nzt+1 )  THEN
    1003                             w_int = w_int_l
    1004                          ELSE
    1005                             w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
    1006                                                                w(kk+1,jj,ii+1) &
    1007                                       + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
    1008                                                              w(kk+1,jj+1,ii+1) &
    1009                                       ) / ( 3.0 * gg )
    1010                             w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
    1011                                               * ( w_int_u - w_int_l )
    1012                          ENDIF
    1013 
    1014 !
    1015 !--                      Change in radius due to collision
    1016                          delta_r = effective_coll_efficiency / 3.0             &
    1017                                    * pi * sl_r3 * ddx * ddy / dz               &
    1018                                    * SQRT( ( u_int - particles(n)%speed_x )**2 &
    1019                                          + ( v_int - particles(n)%speed_y )**2 &
    1020                                          + ( w_int - particles(n)%speed_z )**2 &
    1021                                          ) * dt_3d
    1022 !
    1023 !--                      Change in volume due to collision
    1024                          delta_v = particles(n)%weight_factor                  &
    1025                                    * ( ( particles(n)%radius + delta_r )**3    &
    1026                                        - particles(n)%radius**3 )
    1027 
    1028 !
    1029 !--                      Check if collected particles provide enough LWC for
    1030 !--                      volume change of collector particle
    1031                          IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0 )  THEN
    1032 
    1033                             delta_r = ( ( sl_r3/particles(n)%weight_factor )   &
    1034                                         + particles(n)%radius**3 )**( 1./3. )  &
    1035                                         - particles(n)%radius
    1036 
    1037                             DO  is = n-1, psi, -1
    1038                                IF ( particles(is)%radius < &
    1039                                     particles(n)%radius )  THEN
    1040                                   particles(is)%weight_factor = 0.0
    1041                                   particle_mask(is) = .FALSE.
    1042                                   deleted_particles = deleted_particles + 1
    1043                                ENDIF
    1044                             ENDDO
    1045 
    1046                          ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0 )  THEN
    1047 
    1048                             DO  is = n-1, psi, -1
    1049                                IF ( particles(is)%radius < particles(n)%radius &
    1050                                     .AND.  sl_r3 > 0.0 )  THEN
    1051                                   particles(is)%weight_factor =                &
    1052                                                ( ( particles(is)%weight_factor &
    1053                                                * ( particles(is)%radius**3 ) ) &
    1054                                                - ( delta_v                     &
    1055                                                * particles(is)%weight_factor   &
    1056                                                * ( particles(is)%radius**3 )   &
    1057                                                / sl_r3 ) )                     &
    1058                                                / ( particles(is)%radius**3 )
    1059 
    1060                                   IF ( particles(is)%weight_factor < 0.0 )  THEN
    1061                                      WRITE( message_string, * ) 'negative ',   &
    1062                                                      'weighting factor: ',     &
    1063                                                      particles(is)%weight_factor
    1064                                      CALL message( 'advec_particles', '', 2,   &
    1065                                                    2, -1, 6, 1 )
    1066                                   ENDIF
    1067                                ENDIF
    1068                             ENDDO
    1069                          ENDIF
    1070 
    1071                          particles(n)%radius = particles(n)%radius + delta_r
    1072                          ql_vp(k,j,i) = ql_vp(k,j,i) +               &
    1073                                         particles(n)%weight_factor * &
    1074                                         ( particles(n)%radius**3 )
    1075                       ENDDO
    1076 
    1077                    ENDIF  ! collision kernel
    1078 
    1079                    ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
    1080                                                  * particles(psi)%radius**3
    1081 
    1082 
    1083                 ELSE IF ( prt_count(k,j,i) == 1 )  THEN
    1084 
    1085                    psi = prt_start_index(k,j,i)
    1086                    ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
    1087                                    particles(psi)%radius**3
    1088                 ENDIF
    1089 
    1090 !
    1091 !--             Check if condensation of LWC was conserved during collision
    1092 !--             process
    1093                 IF ( ql_v(k,j,i) /= 0.0 )  THEN
    1094                    IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001  .OR.             &
    1095                         ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 )  THEN
    1096                       WRITE( message_string, * ) 'LWC is not conserved during',&
    1097                                                  ' collision! ',               &
    1098                                                  'LWC after condensation: ',   &
    1099                                                  ql_v(k,j,i),                  &
    1100                                                  ' LWC after collision: ',     &
    1101                                                  ql_vp(k,j,i)
    1102                       CALL message( 'advec_particles', '', 2, 2, -1, 6,  &
    1103                                     1 )
    1104                    ENDIF
    1105                 ENDIF
    1106 
    1107              ENDDO
    1108           ENDDO
    1109        ENDDO
    1110 
    1111        ENDIF   ! collision handling
    1112 
    1113        CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' )
    1114 
    1115     ENDIF   ! cloud droplet handling
    1116 
    1117 
    1118 !
    1119 !-- Particle advection.
    1120 !-- In case of including the SGS velocities, the LES timestep has probably
    1121 !-- to be split into several smaller timesteps because of the Lagrangian
    1122 !-- timescale condition. Because the number of timesteps to be carried out is
    1123 !-- not known at the beginning, these steps are carried out in an infinite loop
    1124 !-- with exit condition.
    1125 !
    1126 !-- If SGS velocities are used, gradients of the TKE have to be calculated and
    1127 !-- boundary conditions have to be set first. Also, horizontally averaged
    1128 !-- profiles of the SGS TKE and the resolved-scale velocity variances are
    1129 !-- needed.
    1130     IF ( use_sgs_for_particles )  THEN
    1131 
    1132 !
    1133 !--    TKE gradient along x and y
    1134        DO  i = nxl, nxr
    1135           DO  j = nys, nyn
    1136              DO  k = nzb, nzt+1
    1137 
    1138                 IF ( k <= nzb_s_inner(j,i-1)  .AND.  &
    1139                      k  > nzb_s_inner(j,i)    .AND.  &
    1140                      k  > nzb_s_inner(j,i+1) )  THEN
    1141                    de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
    1142                                   ( e(k,j,i+1) - e(k,j,i) ) * ddx
    1143                 ELSEIF ( k  > nzb_s_inner(j,i-1)  .AND.  &
    1144                          k  > nzb_s_inner(j,i)    .AND.  &
    1145                          k <= nzb_s_inner(j,i+1) )  THEN
    1146                    de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
    1147                                   ( e(k,j,i) - e(k,j,i-1) ) * ddx
    1148                 ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j,i+1) ) &
    1149                 THEN
    1150                    de_dx(k,j,i) = 0.0
    1151                 ELSEIF ( k < nzb_s_inner(j,i-1)  .AND.  k < nzb_s_inner(j,i) ) &
    1152                 THEN
    1153                    de_dx(k,j,i) = 0.0
    1154                 ELSE
    1155                    de_dx(k,j,i) = sgs_wfu_part * &
    1156                                   ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
    1157                 ENDIF
    1158 
    1159                 IF ( k <= nzb_s_inner(j-1,i)  .AND.  &
    1160                      k  > nzb_s_inner(j,i)    .AND.  &
    1161                      k  > nzb_s_inner(j+1,i) )  THEN
    1162                    de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
    1163                                   ( e(k,j+1,i) - e(k,j,i) ) * ddy
    1164                 ELSEIF ( k  > nzb_s_inner(j-1,i)  .AND.  &
    1165                          k  > nzb_s_inner(j,i)    .AND.  &
    1166                          k <= nzb_s_inner(j+1,i) )  THEN
    1167                    de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
    1168                                   ( e(k,j,i) - e(k,j-1,i) ) * ddy
    1169                 ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j+1,i) ) &
    1170                 THEN
    1171                    de_dy(k,j,i) = 0.0
    1172                 ELSEIF ( k < nzb_s_inner(j-1,i)  .AND.  k < nzb_s_inner(j,i) ) &
    1173                 THEN
    1174                    de_dy(k,j,i) = 0.0
    1175                 ELSE
    1176                    de_dy(k,j,i) = sgs_wfv_part * &
    1177                                   ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
    1178                 ENDIF
    1179 
    1180              ENDDO
    1181           ENDDO
    1182        ENDDO
    1183 
    1184 !
    1185 !--    TKE gradient along z, including bottom and top boundary conditions
    1186        DO  i = nxl, nxr
    1187           DO  j = nys, nyn
    1188 
    1189              DO  k = nzb_s_inner(j,i)+2, nzt-1
    1190                 de_dz(k,j,i) = 2.0 * sgs_wfw_part * &
    1191                                ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
    1192              ENDDO
    1193 
    1194              k = nzb_s_inner(j,i)
    1195              de_dz(nzb:k,j,i)   = 0.0
    1196              de_dz(k+1,j,i) = 2.0 * sgs_wfw_part * ( e(k+2,j,i) - e(k+1,j,i) ) &
    1197                                                  / ( zu(k+2) - zu(k+1) )
    1198              de_dz(nzt,j,i)   = 0.0
    1199              de_dz(nzt+1,j,i) = 0.0
    1200           ENDDO
    1201        ENDDO
    1202 
    1203 !
    1204 !--    Lateral boundary conditions
    1205        CALL exchange_horiz( de_dx, nbgp )
    1206        CALL exchange_horiz( de_dy, nbgp )
    1207        CALL exchange_horiz( de_dz, nbgp )
    1208        CALL exchange_horiz( diss, nbgp  )
    1209 
    1210 !
    1211 !--    Calculate the horizontally averaged profiles of SGS TKE and resolved
    1212 !--    velocity variances (they may have been already calculated in routine
    1213 !--    flow_statistics).
    1214        IF ( .NOT. flow_statistics_called )  THEN
    1215 !
    1216 !--       First calculate horizontally averaged profiles of the horizontal
    1217 !--       velocities.
    1218           sums_l(:,1,0) = 0.0
    1219           sums_l(:,2,0) = 0.0
    1220 
    1221           DO  i = nxl, nxr
    1222              DO  j =  nys, nyn
    1223                 DO  k = nzb_s_outer(j,i), nzt+1
    1224                    sums_l(k,1,0)  = sums_l(k,1,0)  + u(k,j,i)
    1225                    sums_l(k,2,0)  = sums_l(k,2,0)  + v(k,j,i)
    1226                 ENDDO
    1227              ENDDO
    1228           ENDDO
    1229 
    1230 #if defined( __parallel )
    1231 !
    1232 !--       Compute total sum from local sums
    1233           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1234           CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
    1235                               MPI_REAL, MPI_SUM, comm2d, ierr )
    1236           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1237           CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
    1238                               MPI_REAL, MPI_SUM, comm2d, ierr )
    1239 #else
    1240                    sums(:,1) = sums_l(:,1,0)
    1241                    sums(:,2) = sums_l(:,2,0)
    1242 #endif
    1243 
    1244 !
    1245 !--       Final values are obtained by division by the total number of grid
    1246 !--       points used for the summation.
    1247           hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0)   ! u
    1248           hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0)   ! v
    1249 
    1250 !
    1251 !--       Now calculate the profiles of SGS TKE and the resolved-scale
    1252 !--       velocity variances
    1253           sums_l(:,8,0)  = 0.0
    1254           sums_l(:,30,0) = 0.0
    1255           sums_l(:,31,0) = 0.0
    1256           sums_l(:,32,0) = 0.0
    1257           DO  i = nxl, nxr
    1258              DO  j = nys, nyn
    1259                 DO  k = nzb_s_outer(j,i), nzt+1
    1260                    sums_l(k,8,0)  = sums_l(k,8,0)  + e(k,j,i)
    1261                    sums_l(k,30,0) = sums_l(k,30,0) + &
    1262                                     ( u(k,j,i) - hom(k,1,1,0) )**2
    1263                    sums_l(k,31,0) = sums_l(k,31,0) + &
    1264                                     ( v(k,j,i) - hom(k,1,2,0) )**2
    1265                    sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2
    1266                 ENDDO
    1267              ENDDO
    1268           ENDDO
    1269 
    1270 #if defined( __parallel )
    1271 !
    1272 !--       Compute total sum from local sums
    1273           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1274           CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
    1275                               MPI_REAL, MPI_SUM, comm2d, ierr )
    1276           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1277           CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
    1278                               MPI_REAL, MPI_SUM, comm2d, ierr )
    1279           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1280           CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
    1281                               MPI_REAL, MPI_SUM, comm2d, ierr )
    1282           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1283           CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
    1284                               MPI_REAL, MPI_SUM, comm2d, ierr )
    1285 
    1286 #else
    1287           sums(:,8)  = sums_l(:,8,0)
    1288           sums(:,30) = sums_l(:,30,0)
    1289           sums(:,31) = sums_l(:,31,0)
    1290           sums(:,32) = sums_l(:,32,0)
    1291 #endif
    1292 
    1293 !
    1294 !--       Final values are obtained by division by the total number of grid
    1295 !--       points used for the summation.
    1296           hom(:,1,8,0)  = sums(:,8)  / ngp_2dh_outer(:,0)   ! e
    1297           hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0)   ! u*2
    1298           hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0)   ! v*2
    1299           hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0)   ! w*2
    1300 
    1301        ENDIF
    1302 
    1303     ENDIF
     81    ENDIF
     82
     83
     84!
     85!-- Initialize arrays for marking those particles/tails to be deleted after the
     86!-- (sub-) timestep
     87    particle_mask     = .TRUE.
     88    deleted_particles = 0
     89
     90    IF ( use_particle_tails )  THEN
     91       tail_mask = .TRUE.
     92    ENDIF
     93    deleted_tails = 0
     94
    130495
    130596!
     
    1317108    trnp_count_recv_sum = 0
    1318109
     110
     111!
     112!-- Calculate exponential term used in case of particle inertia for each
     113!-- of the particle groups
     114    DO  m = 1, number_of_particle_groups
     115       IF ( particle_groups(m)%density_ratio /= 0.0 )  THEN
     116          particle_groups(m)%exp_arg  =                                        &
     117                    4.5 * particle_groups(m)%density_ratio *                   &
     118                    molecular_viscosity / ( particle_groups(m)%radius )**2
     119          particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * &
     120                                             dt_3d )
     121       ENDIF
     122    ENDDO
     123
     124
     125!
     126!-- Particle (droplet) growth by condensation/evaporation and collision
     127    IF ( cloud_droplets )  THEN
     128
     129!
     130!--    Reset summation arrays
     131       ql_c = 0.0;  ql_v = 0.0;  ql_vp = 0.0
     132
     133!
     134!--    Droplet growth by condensation / evaporation
     135       CALL lpm_droplet_condensation
     136
     137!
     138!--    Particle growth by collision
     139       IF ( collision_kernel /= 'none' )  CALL lpm_droplet_collision
     140
     141    ENDIF
     142
     143
     144!
     145!-- If particle advection includes SGS velocity components, calculate the
     146!-- required SGS quantities (i.e. gradients of the TKE, as well as horizontally
     147!-- averaged profiles of the SGS TKE and the resolved-scale velocity variances)
     148    IF ( use_sgs_for_particles )  CALL lpm_init_sgs_tke
     149
     150
    1319151!
    1320152!-- Initialize the variable storing the total time that a particle has advanced
     
    1322154    particles(1:number_of_particles)%dt_sum = 0.0
    1323155
    1324 !
    1325 !-- Timestep loop.
     156
     157!
     158!-- Timestep loop for particle advection.
    1326159!-- This loop has to be repeated until the advection time of every particle
    1327 !-- (in the total domain!) has reached the LES timestep (dt_3d)
     160!-- (within the total domain!) has reached the LES timestep (dt_3d).
     161!-- In case of including the SGS velocities, the particle timestep may be
     162!-- smaller than the LES timestep (because of the Lagrangian timescale restric-
     163!-- tion) and particles may require to undergo several particle timesteps,
     164!-- before the LES timestep is reached. Because the number of these particle
     165!-- timesteps to be carried out is unknown at first, these steps are carried
     166!-- out in the following infinite loop with exit condition.
    1328167    DO
    1329168
    1330        CALL cpu_log( log_point_s(44), 'advec_part_advec', 'start' )
     169       CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' )
    1331170
    1332171!
     
    1337176       dt_3d_reached_l = .TRUE.
    1338177
    1339        DO  n = 1, number_of_particles
    1340 !
    1341 !--       Move particles only if the LES timestep has not (approximately) been
    1342 !--       reached
    1343           IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 )  CYCLE
    1344 
    1345 !
    1346 !--       Interpolate u velocity-component, determine left, front, bottom
    1347 !--       index of u-array
    1348           i = ( particles(n)%x + 0.5 * dx ) * ddx
    1349           j =   particles(n)%y * ddy
    1350           k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    1351               + offset_ocean_nzt                     ! only exact if equidistant
    1352 
    1353 !
    1354 !--       Interpolation of the velocity components in the xy-plane
    1355           x  = particles(n)%x + ( 0.5 - i ) * dx
    1356           y  = particles(n)%y - j * dy
    1357           aa = x**2          + y**2
    1358           bb = ( dx - x )**2 + y**2
    1359           cc = x**2          + ( dy - y )**2
    1360           dd = ( dx - x )**2 + ( dy - y )**2
    1361           gg = aa + bb + cc + dd
    1362 
    1363           u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
    1364                     + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
    1365                     ) / ( 3.0 * gg ) - u_gtrans
    1366           IF ( k+1 == nzt+1 )  THEN
    1367              u_int = u_int_l
    1368           ELSE
    1369              u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1) &
    1370                     + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1)  &
    1371                     ) / ( 3.0 * gg ) - u_gtrans
    1372              u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1373                                ( u_int_u - u_int_l )
    1374           ENDIF
    1375 
    1376 !
    1377 !--       Same procedure for interpolation of the v velocity-component (adopt
    1378 !--       index k from u velocity-component)
    1379           i =   particles(n)%x * ddx
    1380           j = ( particles(n)%y + 0.5 * dy ) * ddy
    1381 
    1382           x  = particles(n)%x - i * dx
    1383           y  = particles(n)%y + ( 0.5 - j ) * dy
    1384           aa = x**2          + y**2
    1385           bb = ( dx - x )**2 + y**2
    1386           cc = x**2          + ( dy - y )**2
    1387           dd = ( dx - x )**2 + ( dy - y )**2
    1388           gg = aa + bb + cc + dd
    1389 
    1390           v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
    1391                     + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
    1392                     ) / ( 3.0 * gg ) - v_gtrans
    1393           IF ( k+1 == nzt+1 )  THEN
    1394              v_int = v_int_l
    1395           ELSE
    1396              v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1) &
    1397                     + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1)  &
    1398                     ) / ( 3.0 * gg ) - v_gtrans
    1399              v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1400                                ( v_int_u - v_int_l )
    1401           ENDIF
    1402 
    1403 !
    1404 !--       Same procedure for interpolation of the w velocity-component (adopt
    1405 !--       index i from v velocity-component)
    1406           IF ( vertical_particle_advection(particles(n)%group) )  THEN
    1407              j = particles(n)%y * ddy
    1408              k = particles(n)%z / dz + offset_ocean_nzt_m1
    1409 
    1410              x  = particles(n)%x - i * dx
    1411              y  = particles(n)%y - j * dy
    1412              aa = x**2          + y**2
    1413              bb = ( dx - x )**2 + y**2
    1414              cc = x**2          + ( dy - y )**2
    1415              dd = ( dx - x )**2 + ( dy - y )**2
    1416              gg = aa + bb + cc + dd
    1417 
    1418              w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1) &
    1419                     + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1)  &
    1420                        ) / ( 3.0 * gg )
    1421              IF ( k+1 == nzt+1 )  THEN
    1422                 w_int = w_int_l
    1423              ELSE
    1424                 w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
    1425                             ( gg-bb ) * w(k+1,j,i+1) + &
    1426                             ( gg-cc ) * w(k+1,j+1,i) + &
    1427                             ( gg-dd ) * w(k+1,j+1,i+1) &
    1428                            ) / ( 3.0 * gg )
    1429                 w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
    1430                                   ( w_int_u - w_int_l )
    1431              ENDIF
    1432           ELSE
    1433              w_int = 0.0
    1434           ENDIF
    1435 
    1436 !
    1437 !--       Interpolate and calculate quantities needed for calculating the SGS
    1438 !--       velocities
    1439           IF ( use_sgs_for_particles )  THEN
    1440 !
    1441 !--          Interpolate TKE
    1442              i = particles(n)%x * ddx
    1443              j = particles(n)%y * ddy
    1444              k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    1445                  + offset_ocean_nzt                      ! only exact if eq.dist
    1446 
    1447              IF ( topography == 'flat' )  THEN
    1448 
    1449                 x  = particles(n)%x - i * dx
    1450                 y  = particles(n)%y - j * dy
    1451                 aa = x**2          + y**2
    1452                 bb = ( dx - x )**2 + y**2
    1453                 cc = x**2          + ( dy - y )**2
    1454                 dd = ( dx - x )**2 + ( dy - y )**2
    1455                 gg = aa + bb + cc + dd
    1456 
    1457                 e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
    1458                           + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
    1459                           ) / ( 3.0 * gg )
    1460 
    1461                 IF ( k+1 == nzt+1 )  THEN
    1462                    e_int = e_int_l
    1463                 ELSE
    1464                    e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
    1465                                ( gg - bb ) * e(k+1,j,i+1) + &
    1466                                ( gg - cc ) * e(k+1,j+1,i) + &
    1467                                ( gg - dd ) * e(k+1,j+1,i+1) &
    1468                              ) / ( 3.0 * gg )
    1469                    e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1470                                      ( e_int_u - e_int_l )
    1471                 ENDIF
    1472 
    1473 !
    1474 !--             Interpolate the TKE gradient along x (adopt incides i,j,k and
    1475 !--             all position variables from above (TKE))
    1476                 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
    1477                                 ( gg - bb ) * de_dx(k,j,i+1) + &
    1478                                 ( gg - cc ) * de_dx(k,j+1,i) + &
    1479                                 ( gg - dd ) * de_dx(k,j+1,i+1) &
    1480                               ) / ( 3.0 * gg )
    1481 
    1482                 IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1483                    de_dx_int = de_dx_int_l
    1484                 ELSE
    1485                    de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
    1486                                    ( gg - bb ) * de_dx(k+1,j,i+1) + &
    1487                                    ( gg - cc ) * de_dx(k+1,j+1,i) + &
    1488                                    ( gg - dd ) * de_dx(k+1,j+1,i+1) &
    1489                                  ) / ( 3.0 * gg )
    1490                    de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1491                                              ( de_dx_int_u - de_dx_int_l )
    1492                 ENDIF
    1493 
    1494 !
    1495 !--             Interpolate the TKE gradient along y
    1496                 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
    1497                                 ( gg - bb ) * de_dy(k,j,i+1) + &
    1498                                 ( gg - cc ) * de_dy(k,j+1,i) + &
    1499                                 ( gg - dd ) * de_dy(k,j+1,i+1) &
    1500                               ) / ( 3.0 * gg )
    1501                 IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1502                    de_dy_int = de_dy_int_l
    1503                 ELSE
    1504                    de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
    1505                                    ( gg - bb ) * de_dy(k+1,j,i+1) + &
    1506                                    ( gg - cc ) * de_dy(k+1,j+1,i) + &
    1507                                    ( gg - dd ) * de_dy(k+1,j+1,i+1) &
    1508                                  ) / ( 3.0 * gg )
    1509                    de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1510                                              ( de_dy_int_u - de_dy_int_l )
    1511                 ENDIF
    1512 
    1513 !
    1514 !--             Interpolate the TKE gradient along z
    1515                 IF ( particles(n)%z < 0.5 * dz ) THEN
    1516                    de_dz_int = 0.0
    1517                 ELSE
    1518                    de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
    1519                                    ( gg - bb ) * de_dz(k,j,i+1) + &
    1520                                    ( gg - cc ) * de_dz(k,j+1,i) + &
    1521                                    ( gg - dd ) * de_dz(k,j+1,i+1) &
    1522                                  ) / ( 3.0 * gg )
    1523 
    1524                    IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1525                       de_dz_int = de_dz_int_l
    1526                    ELSE
    1527                       de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
    1528                                       ( gg - bb ) * de_dz(k+1,j,i+1) + &
    1529                                       ( gg - cc ) * de_dz(k+1,j+1,i) + &
    1530                                       ( gg - dd ) * de_dz(k+1,j+1,i+1) &
    1531                                     ) / ( 3.0 * gg )
    1532                       de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1533                                                 ( de_dz_int_u - de_dz_int_l )
    1534                    ENDIF
    1535                 ENDIF
    1536 
    1537 !
    1538 !--             Interpolate the dissipation of TKE
    1539                 diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
    1540                                ( gg - bb ) * diss(k,j,i+1) + &
    1541                                ( gg - cc ) * diss(k,j+1,i) + &
    1542                                ( gg - dd ) * diss(k,j+1,i+1) &
    1543                               ) / ( 3.0 * gg )
    1544 
    1545                 IF ( k+1 == nzt+1 )  THEN
    1546                    diss_int = diss_int_l
    1547                 ELSE
    1548                    diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
    1549                                   ( gg - bb ) * diss(k+1,j,i+1) + &
    1550                                   ( gg - cc ) * diss(k+1,j+1,i) + &
    1551                                   ( gg - dd ) * diss(k+1,j+1,i+1) &
    1552                                  ) / ( 3.0 * gg )
    1553                    diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1554                                            ( diss_int_u - diss_int_l )
    1555                 ENDIF
    1556 
    1557              ELSE
    1558 
    1559 !
    1560 !--             In case that there are buildings it has to be determined
    1561 !--             how many of the gridpoints defining the particle box are
    1562 !--             situated within a building
    1563 !--             gp_outside_of_building(1): i,j,k
    1564 !--             gp_outside_of_building(2): i,j+1,k
    1565 !--             gp_outside_of_building(3): i,j,k+1
    1566 !--             gp_outside_of_building(4): i,j+1,k+1
    1567 !--             gp_outside_of_building(5): i+1,j,k
    1568 !--             gp_outside_of_building(6): i+1,j+1,k
    1569 !--             gp_outside_of_building(7): i+1,j,k+1
    1570 !--             gp_outside_of_building(8): i+1,j+1,k+1
    1571 
    1572                 gp_outside_of_building = 0
    1573                 location = 0.0
    1574                 num_gp = 0
    1575 
    1576                 IF ( k > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
    1577                    num_gp = num_gp + 1
    1578                    gp_outside_of_building(1) = 1
    1579                    location(num_gp,1) = i * dx
    1580                    location(num_gp,2) = j * dy
    1581                    location(num_gp,3) = k * dz - 0.5 * dz
    1582                    ei(num_gp)     = e(k,j,i)
    1583                    dissi(num_gp)  = diss(k,j,i)
    1584                    de_dxi(num_gp) = de_dx(k,j,i)
    1585                    de_dyi(num_gp) = de_dy(k,j,i)
    1586                    de_dzi(num_gp) = de_dz(k,j,i)
    1587                 ENDIF
    1588 
    1589                 IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
    1590                 THEN
    1591                    num_gp = num_gp + 1
    1592                    gp_outside_of_building(2) = 1
    1593                    location(num_gp,1) = i * dx
    1594                    location(num_gp,2) = (j+1) * dy
    1595                    location(num_gp,3) = k * dz - 0.5 * dz
    1596                    ei(num_gp)     = e(k,j+1,i)
    1597                    dissi(num_gp)  = diss(k,j+1,i)
    1598                    de_dxi(num_gp) = de_dx(k,j+1,i)
    1599                    de_dyi(num_gp) = de_dy(k,j+1,i)
    1600                    de_dzi(num_gp) = de_dz(k,j+1,i)
    1601                 ENDIF
    1602 
    1603                 IF ( k+1 > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
    1604                    num_gp = num_gp + 1
    1605                    gp_outside_of_building(3) = 1
    1606                    location(num_gp,1) = i * dx
    1607                    location(num_gp,2) = j * dy
    1608                    location(num_gp,3) = (k+1) * dz - 0.5 * dz
    1609                    ei(num_gp)     = e(k+1,j,i)
    1610                    dissi(num_gp)  = diss(k+1,j,i)
    1611                    de_dxi(num_gp) = de_dx(k+1,j,i)
    1612                    de_dyi(num_gp) = de_dy(k+1,j,i)
    1613                    de_dzi(num_gp) = de_dz(k+1,j,i)
    1614                 ENDIF
    1615 
    1616                 IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
    1617                 THEN
    1618                    num_gp = num_gp + 1
    1619                    gp_outside_of_building(4) = 1
    1620                    location(num_gp,1) = i * dx
    1621                    location(num_gp,2) = (j+1) * dy
    1622                    location(num_gp,3) = (k+1) * dz - 0.5 * dz
    1623                    ei(num_gp)     = e(k+1,j+1,i)
    1624                    dissi(num_gp)  = diss(k+1,j+1,i)
    1625                    de_dxi(num_gp) = de_dx(k+1,j+1,i)
    1626                    de_dyi(num_gp) = de_dy(k+1,j+1,i)
    1627                    de_dzi(num_gp) = de_dz(k+1,j+1,i)
    1628                 ENDIF
    1629 
    1630                 IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
    1631                 THEN
    1632                    num_gp = num_gp + 1
    1633                    gp_outside_of_building(5) = 1
    1634                    location(num_gp,1) = (i+1) * dx
    1635                    location(num_gp,2) = j * dy
    1636                    location(num_gp,3) = k * dz - 0.5 * dz
    1637                    ei(num_gp)     = e(k,j,i+1)
    1638                    dissi(num_gp)  = diss(k,j,i+1)
    1639                    de_dxi(num_gp) = de_dx(k,j,i+1)
    1640                    de_dyi(num_gp) = de_dy(k,j,i+1)
    1641                    de_dzi(num_gp) = de_dz(k,j,i+1)
    1642                 ENDIF
    1643 
    1644                 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
    1645                 THEN
    1646                    num_gp = num_gp + 1
    1647                    gp_outside_of_building(6) = 1
    1648                    location(num_gp,1) = (i+1) * dx
    1649                    location(num_gp,2) = (j+1) * dy
    1650                    location(num_gp,3) = k * dz - 0.5 * dz
    1651                    ei(num_gp)     = e(k,j+1,i+1)
    1652                    dissi(num_gp)  = diss(k,j+1,i+1)
    1653                    de_dxi(num_gp) = de_dx(k,j+1,i+1)
    1654                    de_dyi(num_gp) = de_dy(k,j+1,i+1)
    1655                    de_dzi(num_gp) = de_dz(k,j+1,i+1)
    1656                 ENDIF
    1657 
    1658                 IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
    1659                 THEN
    1660                    num_gp = num_gp + 1
    1661                    gp_outside_of_building(7) = 1
    1662                    location(num_gp,1) = (i+1) * dx
    1663                    location(num_gp,2) = j * dy
    1664                    location(num_gp,3) = (k+1) * dz - 0.5 * dz
    1665                    ei(num_gp)     = e(k+1,j,i+1)
    1666                    dissi(num_gp)  = diss(k+1,j,i+1)
    1667                    de_dxi(num_gp) = de_dx(k+1,j,i+1)
    1668                    de_dyi(num_gp) = de_dy(k+1,j,i+1)
    1669                    de_dzi(num_gp) = de_dz(k+1,j,i+1)
    1670                 ENDIF
    1671 
    1672                 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
    1673                 THEN
    1674                    num_gp = num_gp + 1
    1675                    gp_outside_of_building(8) = 1
    1676                    location(num_gp,1) = (i+1) * dx
    1677                    location(num_gp,2) = (j+1) * dy
    1678                    location(num_gp,3) = (k+1) * dz - 0.5 * dz
    1679                    ei(num_gp)     = e(k+1,j+1,i+1)
    1680                    dissi(num_gp)  = diss(k+1,j+1,i+1)
    1681                    de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    1682                    de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    1683                    de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    1684                 ENDIF
    1685 
    1686 !
    1687 !--             If all gridpoints are situated outside of a building, then the
    1688 !--             ordinary interpolation scheme can be used.
    1689                 IF ( num_gp == 8 )  THEN
    1690 
    1691                    x  = particles(n)%x - i * dx
    1692                    y  = particles(n)%y - j * dy
    1693                    aa = x**2          + y**2
    1694                    bb = ( dx - x )**2 + y**2
    1695                    cc = x**2          + ( dy - y )**2
    1696                    dd = ( dx - x )**2 + ( dy - y )**2
    1697                    gg = aa + bb + cc + dd
    1698 
    1699                    e_int_l = (( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1) &
    1700                             + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
    1701                              ) / ( 3.0 * gg )
    1702 
    1703                    IF ( k+1 == nzt+1 )  THEN
    1704                       e_int = e_int_l
    1705                    ELSE
    1706                       e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
    1707                                   ( gg - bb ) * e(k+1,j,i+1) + &
    1708                                   ( gg - cc ) * e(k+1,j+1,i) + &
    1709                                   ( gg - dd ) * e(k+1,j+1,i+1) &
    1710                                 ) / ( 3.0 * gg )
    1711                       e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1712                                         ( e_int_u - e_int_l )
    1713                    ENDIF
    1714 
    1715 !
    1716 !--                Interpolate the TKE gradient along x (adopt incides i,j,k
    1717 !--                and all position variables from above (TKE))
    1718                    de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
    1719                                    ( gg - bb ) * de_dx(k,j,i+1) + &
    1720                                    ( gg - cc ) * de_dx(k,j+1,i) + &
    1721                                    ( gg - dd ) * de_dx(k,j+1,i+1) &
    1722                                  ) / ( 3.0 * gg )
    1723 
    1724                    IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1725                       de_dx_int = de_dx_int_l
    1726                    ELSE
    1727                       de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
    1728                                       ( gg - bb ) * de_dx(k+1,j,i+1) + &
    1729                                       ( gg - cc ) * de_dx(k+1,j+1,i) + &
    1730                                       ( gg - dd ) * de_dx(k+1,j+1,i+1) &
    1731                                     ) / ( 3.0 * gg )
    1732                       de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
    1733                                               dz * ( de_dx_int_u - de_dx_int_l )
    1734                    ENDIF
    1735 
    1736 !
    1737 !--                Interpolate the TKE gradient along y
    1738                    de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
    1739                                    ( gg - bb ) * de_dy(k,j,i+1) + &
    1740                                    ( gg - cc ) * de_dy(k,j+1,i) + &
    1741                                    ( gg - dd ) * de_dy(k,j+1,i+1) &
    1742                                  ) / ( 3.0 * gg )
    1743                    IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1744                       de_dy_int = de_dy_int_l
    1745                    ELSE
    1746                       de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
    1747                                       ( gg - bb ) * de_dy(k+1,j,i+1) + &
    1748                                       ( gg - cc ) * de_dy(k+1,j+1,i) + &
    1749                                       ( gg - dd ) * de_dy(k+1,j+1,i+1) &
    1750                                     ) / ( 3.0 * gg )
    1751                       de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
    1752                                               dz * ( de_dy_int_u - de_dy_int_l )
    1753                    ENDIF
    1754 
    1755 !
    1756 !--                Interpolate the TKE gradient along z
    1757                    IF ( particles(n)%z < 0.5 * dz )  THEN
    1758                       de_dz_int = 0.0
    1759                    ELSE
    1760                       de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
    1761                                       ( gg - bb ) * de_dz(k,j,i+1) + &
    1762                                       ( gg - cc ) * de_dz(k,j+1,i) + &
    1763                                       ( gg - dd ) * de_dz(k,j+1,i+1) &
    1764                                     ) / ( 3.0 * gg )
    1765 
    1766                       IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1767                          de_dz_int = de_dz_int_l
    1768                       ELSE
    1769                          de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
    1770                                          ( gg - bb ) * de_dz(k+1,j,i+1) + &
    1771                                          ( gg - cc ) * de_dz(k+1,j+1,i) + &
    1772                                          ( gg - dd ) * de_dz(k+1,j+1,i+1) &
    1773                                        ) / ( 3.0 * gg )
    1774                          de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
    1775                                               dz * ( de_dz_int_u - de_dz_int_l )
    1776                       ENDIF
    1777                    ENDIF
    1778 
    1779 !
    1780 !--                Interpolate the dissipation of TKE
    1781                    diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
    1782                                   ( gg - bb ) * diss(k,j,i+1) + &
    1783                                   ( gg - cc ) * diss(k,j+1,i) + &
    1784                                   ( gg - dd ) * diss(k,j+1,i+1) &
    1785                                 ) / ( 3.0 * gg )
    1786 
    1787                    IF ( k+1 == nzt+1 )  THEN
    1788                       diss_int = diss_int_l
    1789                    ELSE
    1790                       diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
    1791                                      ( gg - bb ) * diss(k+1,j,i+1) + &
    1792                                      ( gg - cc ) * diss(k+1,j+1,i) + &
    1793                                      ( gg - dd ) * diss(k+1,j+1,i+1) &
    1794                                    ) / ( 3.0 * gg )
    1795                       diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
    1796                                               ( diss_int_u - diss_int_l )
    1797                    ENDIF
    1798 
    1799                 ELSE
    1800 
    1801 !
    1802 !--                If wall between gridpoint 1 and gridpoint 5, then
    1803 !--                Neumann boundary condition has to be applied
    1804                    IF ( gp_outside_of_building(1) == 1  .AND. &
    1805                         gp_outside_of_building(5) == 0 )  THEN
    1806                       num_gp = num_gp + 1
    1807                       location(num_gp,1) = i * dx + 0.5 * dx
    1808                       location(num_gp,2) = j * dy
    1809                       location(num_gp,3) = k * dz - 0.5 * dz
    1810                       ei(num_gp)     = e(k,j,i)
    1811                       dissi(num_gp)  = diss(k,j,i)
    1812                       de_dxi(num_gp) = 0.0
    1813                       de_dyi(num_gp) = de_dy(k,j,i)
    1814                       de_dzi(num_gp) = de_dz(k,j,i)
    1815                    ENDIF
    1816 
    1817                    IF ( gp_outside_of_building(5) == 1  .AND. &
    1818                       gp_outside_of_building(1) == 0 )  THEN
    1819                       num_gp = num_gp + 1
    1820                       location(num_gp,1) = i * dx + 0.5 * dx
    1821                       location(num_gp,2) = j * dy
    1822                       location(num_gp,3) = k * dz - 0.5 * dz
    1823                       ei(num_gp)     = e(k,j,i+1)
    1824                       dissi(num_gp)  = diss(k,j,i+1)
    1825                       de_dxi(num_gp) = 0.0
    1826                       de_dyi(num_gp) = de_dy(k,j,i+1)
    1827                       de_dzi(num_gp) = de_dz(k,j,i+1)
    1828                    ENDIF
    1829 
    1830 !
    1831 !--                If wall between gridpoint 5 and gridpoint 6, then
    1832 !--                then Neumann boundary condition has to be applied
    1833                    IF ( gp_outside_of_building(5) == 1  .AND. &
    1834                         gp_outside_of_building(6) == 0 )  THEN
    1835                       num_gp = num_gp + 1
    1836                       location(num_gp,1) = (i+1) * dx
    1837                       location(num_gp,2) = j * dy + 0.5 * dy
    1838                       location(num_gp,3) = k * dz - 0.5 * dz
    1839                       ei(num_gp)     = e(k,j,i+1)
    1840                       dissi(num_gp)  = diss(k,j,i+1)
    1841                       de_dxi(num_gp) = de_dx(k,j,i+1)
    1842                       de_dyi(num_gp) = 0.0
    1843                       de_dzi(num_gp) = de_dz(k,j,i+1)
    1844                    ENDIF
    1845 
    1846                    IF ( gp_outside_of_building(6) == 1  .AND. &
    1847                         gp_outside_of_building(5) == 0 )  THEN
    1848                       num_gp = num_gp + 1
    1849                       location(num_gp,1) = (i+1) * dx
    1850                       location(num_gp,2) = j * dy + 0.5 * dy
    1851                       location(num_gp,3) = k * dz - 0.5 * dz
    1852                       ei(num_gp)     = e(k,j+1,i+1)
    1853                       dissi(num_gp)  = diss(k,j+1,i+1)
    1854                       de_dxi(num_gp) = de_dx(k,j+1,i+1)
    1855                       de_dyi(num_gp) = 0.0
    1856                       de_dzi(num_gp) = de_dz(k,j+1,i+1)
    1857                    ENDIF
    1858 
    1859 !
    1860 !--                If wall between gridpoint 2 and gridpoint 6, then
    1861 !--                Neumann boundary condition has to be applied
    1862                    IF ( gp_outside_of_building(2) == 1  .AND. &
    1863                         gp_outside_of_building(6) == 0 )  THEN
    1864                       num_gp = num_gp + 1
    1865                       location(num_gp,1) = i * dx + 0.5 * dx
    1866                       location(num_gp,2) = (j+1) * dy
    1867                       location(num_gp,3) = k * dz - 0.5 * dz
    1868                       ei(num_gp)     = e(k,j+1,i)
    1869                       dissi(num_gp)  = diss(k,j+1,i)
    1870                       de_dxi(num_gp) = 0.0
    1871                       de_dyi(num_gp) = de_dy(k,j+1,i)
    1872                       de_dzi(num_gp) = de_dz(k,j+1,i)
    1873                    ENDIF
    1874 
    1875                    IF ( gp_outside_of_building(6) == 1  .AND. &
    1876                         gp_outside_of_building(2) == 0 )  THEN
    1877                       num_gp = num_gp + 1
    1878                       location(num_gp,1) = i * dx + 0.5 * dx
    1879                       location(num_gp,2) = (j+1) * dy
    1880                       location(num_gp,3) = k * dz - 0.5 * dz
    1881                       ei(num_gp)     = e(k,j+1,i+1)
    1882                       dissi(num_gp)  = diss(k,j+1,i+1)
    1883                       de_dxi(num_gp) = 0.0
    1884                       de_dyi(num_gp) = de_dy(k,j+1,i+1)
    1885                       de_dzi(num_gp) = de_dz(k,j+1,i+1)
    1886                    ENDIF
    1887 
    1888 !
    1889 !--                If wall between gridpoint 1 and gridpoint 2, then
    1890 !--                Neumann boundary condition has to be applied
    1891                    IF ( gp_outside_of_building(1) == 1  .AND. &
    1892                         gp_outside_of_building(2) == 0 )  THEN
    1893                       num_gp = num_gp + 1
    1894                       location(num_gp,1) = i * dx
    1895                       location(num_gp,2) = j * dy + 0.5 * dy
    1896                       location(num_gp,3) = k * dz - 0.5 * dz
    1897                       ei(num_gp)     = e(k,j,i)
    1898                       dissi(num_gp)  = diss(k,j,i)
    1899                       de_dxi(num_gp) = de_dx(k,j,i)
    1900                       de_dyi(num_gp) = 0.0
    1901                       de_dzi(num_gp) = de_dz(k,j,i)
    1902                    ENDIF
    1903 
    1904                    IF ( gp_outside_of_building(2) == 1  .AND. &
    1905                         gp_outside_of_building(1) == 0 )  THEN
    1906                       num_gp = num_gp + 1
    1907                       location(num_gp,1) = i * dx
    1908                       location(num_gp,2) = j * dy + 0.5 * dy
    1909                       location(num_gp,3) = k * dz - 0.5 * dz
    1910                       ei(num_gp)     = e(k,j+1,i)
    1911                       dissi(num_gp)  = diss(k,j+1,i)
    1912                       de_dxi(num_gp) = de_dx(k,j+1,i)
    1913                       de_dyi(num_gp) = 0.0
    1914                       de_dzi(num_gp) = de_dz(k,j+1,i)
    1915                    ENDIF
    1916 
    1917 !
    1918 !--                If wall between gridpoint 3 and gridpoint 7, then
    1919 !--                Neumann boundary condition has to be applied
    1920                    IF ( gp_outside_of_building(3) == 1  .AND. &
    1921                         gp_outside_of_building(7) == 0 )  THEN
    1922                       num_gp = num_gp + 1
    1923                       location(num_gp,1) = i * dx + 0.5 * dx
    1924                       location(num_gp,2) = j * dy
    1925                       location(num_gp,3) = k * dz + 0.5 * dz
    1926                       ei(num_gp)     = e(k+1,j,i)
    1927                       dissi(num_gp)  = diss(k+1,j,i)
    1928                       de_dxi(num_gp) = 0.0
    1929                       de_dyi(num_gp) = de_dy(k+1,j,i)
    1930                       de_dzi(num_gp) = de_dz(k+1,j,i)
    1931                    ENDIF
    1932 
    1933                    IF ( gp_outside_of_building(7) == 1  .AND. &
    1934                         gp_outside_of_building(3) == 0 )  THEN
    1935                       num_gp = num_gp + 1
    1936                       location(num_gp,1) = i * dx + 0.5 * dx
    1937                       location(num_gp,2) = j * dy
    1938                       location(num_gp,3) = k * dz + 0.5 * dz
    1939                       ei(num_gp)     = e(k+1,j,i+1)
    1940                       dissi(num_gp)  = diss(k+1,j,i+1)
    1941                       de_dxi(num_gp) = 0.0
    1942                       de_dyi(num_gp) = de_dy(k+1,j,i+1)
    1943                       de_dzi(num_gp) = de_dz(k+1,j,i+1)
    1944                    ENDIF
    1945 
    1946 !
    1947 !--                If wall between gridpoint 7 and gridpoint 8, then
    1948 !--                Neumann boundary condition has to be applied
    1949                    IF ( gp_outside_of_building(7) == 1  .AND. &
    1950                         gp_outside_of_building(8) == 0 )  THEN
    1951                       num_gp = num_gp + 1
    1952                       location(num_gp,1) = (i+1) * dx
    1953                       location(num_gp,2) = j * dy + 0.5 * dy
    1954                       location(num_gp,3) = k * dz + 0.5 * dz
    1955                       ei(num_gp)     = e(k+1,j,i+1)
    1956                       dissi(num_gp)  = diss(k+1,j,i+1)
    1957                       de_dxi(num_gp) = de_dx(k+1,j,i+1)
    1958                       de_dyi(num_gp) = 0.0
    1959                       de_dzi(num_gp) = de_dz(k+1,j,i+1)
    1960                    ENDIF
    1961 
    1962                    IF ( gp_outside_of_building(8) == 1  .AND. &
    1963                         gp_outside_of_building(7) == 0 )  THEN
    1964                       num_gp = num_gp + 1
    1965                       location(num_gp,1) = (i+1) * dx
    1966                       location(num_gp,2) = j * dy + 0.5 * dy
    1967                       location(num_gp,3) = k * dz + 0.5 * dz
    1968                       ei(num_gp)     = e(k+1,j+1,i+1)
    1969                       dissi(num_gp)  = diss(k+1,j+1,i+1)
    1970                       de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    1971                       de_dyi(num_gp) = 0.0
    1972                       de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    1973                    ENDIF
    1974 
    1975 !
    1976 !--                If wall between gridpoint 4 and gridpoint 8, then
    1977 !--                Neumann boundary condition has to be applied
    1978                    IF ( gp_outside_of_building(4) == 1  .AND. &
    1979                         gp_outside_of_building(8) == 0 )  THEN
    1980                       num_gp = num_gp + 1
    1981                       location(num_gp,1) = i * dx + 0.5 * dx
    1982                       location(num_gp,2) = (j+1) * dy
    1983                       location(num_gp,3) = k * dz + 0.5 * dz
    1984                       ei(num_gp)     = e(k+1,j+1,i)
    1985                       dissi(num_gp)  = diss(k+1,j+1,i)
    1986                       de_dxi(num_gp) = 0.0
    1987                       de_dyi(num_gp) = de_dy(k+1,j+1,i)
    1988                       de_dzi(num_gp) = de_dz(k+1,j+1,i)
    1989                    ENDIF
    1990 
    1991                    IF ( gp_outside_of_building(8) == 1  .AND. &
    1992                         gp_outside_of_building(4) == 0 )  THEN
    1993                       num_gp = num_gp + 1
    1994                       location(num_gp,1) = i * dx + 0.5 * dx
    1995                       location(num_gp,2) = (j+1) * dy
    1996                       location(num_gp,3) = k * dz + 0.5 * dz
    1997                       ei(num_gp)     = e(k+1,j+1,i+1)
    1998                       dissi(num_gp)  = diss(k+1,j+1,i+1)
    1999                       de_dxi(num_gp) = 0.0
    2000                       de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    2001                       de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    2002                    ENDIF
    2003 
    2004 !
    2005 !--                If wall between gridpoint 3 and gridpoint 4, then
    2006 !--                Neumann boundary condition has to be applied
    2007                    IF ( gp_outside_of_building(3) == 1  .AND. &
    2008                         gp_outside_of_building(4) == 0 )  THEN
    2009                       num_gp = num_gp + 1
    2010                       location(num_gp,1) = i * dx
    2011                       location(num_gp,2) = j * dy + 0.5 * dy
    2012                       location(num_gp,3) = k * dz + 0.5 * dz
    2013                       ei(num_gp)     = e(k+1,j,i)
    2014                       dissi(num_gp)  = diss(k+1,j,i)
    2015                       de_dxi(num_gp) = de_dx(k+1,j,i)
    2016                       de_dyi(num_gp) = 0.0
    2017                       de_dzi(num_gp) = de_dz(k+1,j,i)
    2018                    ENDIF
    2019 
    2020                    IF ( gp_outside_of_building(4) == 1  .AND. &
    2021                         gp_outside_of_building(3) == 0 )  THEN
    2022                       num_gp = num_gp + 1
    2023                       location(num_gp,1) = i * dx
    2024                       location(num_gp,2) = j * dy + 0.5 * dy
    2025                       location(num_gp,3) = k * dz + 0.5 * dz
    2026                       ei(num_gp)     = e(k+1,j+1,i)
    2027                       dissi(num_gp)  = diss(k+1,j+1,i)
    2028                       de_dxi(num_gp) = de_dx(k+1,j+1,i)
    2029                       de_dyi(num_gp) = 0.0
    2030                       de_dzi(num_gp) = de_dz(k+1,j+1,i)
    2031                    ENDIF
    2032 
    2033 !
    2034 !--                If wall between gridpoint 1 and gridpoint 3, then
    2035 !--                Neumann boundary condition has to be applied
    2036 !--                (only one case as only building beneath is possible)
    2037                    IF ( gp_outside_of_building(1) == 0  .AND. &
    2038                         gp_outside_of_building(3) == 1 )  THEN
    2039                       num_gp = num_gp + 1
    2040                       location(num_gp,1) = i * dx
    2041                       location(num_gp,2) = j * dy
    2042                       location(num_gp,3) = k * dz
    2043                       ei(num_gp)     = e(k+1,j,i)
    2044                       dissi(num_gp)  = diss(k+1,j,i)
    2045                       de_dxi(num_gp) = de_dx(k+1,j,i)
    2046                       de_dyi(num_gp) = de_dy(k+1,j,i)
    2047                       de_dzi(num_gp) = 0.0
    2048                    ENDIF
    2049 
    2050 !
    2051 !--                If wall between gridpoint 5 and gridpoint 7, then
    2052 !--                Neumann boundary condition has to be applied
    2053 !--                (only one case as only building beneath is possible)
    2054                    IF ( gp_outside_of_building(5) == 0  .AND. &
    2055                         gp_outside_of_building(7) == 1 )  THEN
    2056                       num_gp = num_gp + 1
    2057                       location(num_gp,1) = (i+1) * dx
    2058                       location(num_gp,2) = j * dy
    2059                       location(num_gp,3) = k * dz
    2060                       ei(num_gp)     = e(k+1,j,i+1)
    2061                       dissi(num_gp)  = diss(k+1,j,i+1)
    2062                       de_dxi(num_gp) = de_dx(k+1,j,i+1)
    2063                       de_dyi(num_gp) = de_dy(k+1,j,i+1)
    2064                       de_dzi(num_gp) = 0.0
    2065                    ENDIF
    2066 
    2067 !
    2068 !--                If wall between gridpoint 2 and gridpoint 4, then
    2069 !--                Neumann boundary condition has to be applied
    2070 !--                (only one case as only building beneath is possible)
    2071                    IF ( gp_outside_of_building(2) == 0  .AND. &
    2072                         gp_outside_of_building(4) == 1 )  THEN
    2073                       num_gp = num_gp + 1
    2074                       location(num_gp,1) = i * dx
    2075                       location(num_gp,2) = (j+1) * dy
    2076                       location(num_gp,3) = k * dz
    2077                       ei(num_gp)     = e(k+1,j+1,i)
    2078                       dissi(num_gp)  = diss(k+1,j+1,i)
    2079                       de_dxi(num_gp) = de_dx(k+1,j+1,i)
    2080                       de_dyi(num_gp) = de_dy(k+1,j+1,i)
    2081                       de_dzi(num_gp) = 0.0
    2082                    ENDIF
    2083 
    2084 !
    2085 !--                If wall between gridpoint 6 and gridpoint 8, then
    2086 !--                Neumann boundary condition has to be applied
    2087 !--                (only one case as only building beneath is possible)
    2088                    IF ( gp_outside_of_building(6) == 0  .AND. &
    2089                         gp_outside_of_building(8) == 1 )  THEN
    2090                       num_gp = num_gp + 1
    2091                       location(num_gp,1) = (i+1) * dx
    2092                       location(num_gp,2) = (j+1) * dy
    2093                       location(num_gp,3) = k * dz
    2094                       ei(num_gp)     = e(k+1,j+1,i+1)
    2095                       dissi(num_gp)  = diss(k+1,j+1,i+1)
    2096                       de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    2097                       de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    2098                       de_dzi(num_gp) = 0.0
    2099                    ENDIF
    2100 
    2101 !
    2102 !--                Carry out the interpolation
    2103                    IF ( num_gp == 1 )  THEN
    2104 !
    2105 !--                   If only one of the gridpoints is situated outside of the
    2106 !--                   building, it follows that the values at the particle
    2107 !--                   location are the same as the gridpoint values
    2108                       e_int     = ei(num_gp)
    2109                       diss_int  = dissi(num_gp)
    2110                       de_dx_int = de_dxi(num_gp)
    2111                       de_dy_int = de_dyi(num_gp)
    2112                       de_dz_int = de_dzi(num_gp)
    2113                    ELSE IF ( num_gp > 1 )  THEN
    2114 
    2115                       d_sum = 0.0
    2116 !
    2117 !--                   Evaluation of the distances between the gridpoints
    2118 !--                   contributing to the interpolated values, and the particle
    2119 !--                   location
    2120                       DO agp = 1, num_gp
    2121                          d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
    2122                                       + ( particles(n)%y-location(agp,2) )**2  &
    2123                                       + ( particles(n)%z-location(agp,3) )**2
    2124                          d_sum        = d_sum + d_gp_pl(agp)
    2125                       ENDDO
    2126 
    2127 !
    2128 !--                   Finally the interpolation can be carried out
    2129                       e_int     = 0.0
    2130                       diss_int  = 0.0
    2131                       de_dx_int = 0.0
    2132                       de_dy_int = 0.0
    2133                       de_dz_int = 0.0
    2134                       DO agp = 1, num_gp
    2135                          e_int     = e_int     + ( d_sum - d_gp_pl(agp) ) * &
    2136                                                 ei(agp) / ( (num_gp-1) * d_sum )
    2137                          diss_int  = diss_int  + ( d_sum - d_gp_pl(agp) ) * &
    2138                                              dissi(agp) / ( (num_gp-1) * d_sum )
    2139                          de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
    2140                                             de_dxi(agp) / ( (num_gp-1) * d_sum )
    2141                          de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
    2142                                             de_dyi(agp) / ( (num_gp-1) * d_sum )
    2143                          de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
    2144                                             de_dzi(agp) / ( (num_gp-1) * d_sum )
    2145                       ENDDO
    2146 
    2147                    ENDIF
    2148 
    2149                 ENDIF
    2150 
    2151              ENDIF 
    2152 
    2153 !
    2154 !--          Vertically interpolate the horizontally averaged SGS TKE and
    2155 !--          resolved-scale velocity variances and use the interpolated values
    2156 !--          to calculate the coefficient fs, which is a measure of the ratio
    2157 !--          of the subgrid-scale turbulent kinetic energy to the total amount
    2158 !--          of turbulent kinetic energy.
    2159              IF ( k == 0 )  THEN
    2160                 e_mean_int = hom(0,1,8,0)
    2161              ELSE
    2162                 e_mean_int = hom(k,1,8,0) +                                    &
    2163                                            ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
    2164                                            ( zu(k+1) - zu(k) ) *               &
    2165                                            ( particles(n)%z - zu(k) )
    2166              ENDIF
    2167 
    2168              kw = particles(n)%z / dz
    2169 
    2170              IF ( k == 0 )  THEN
    2171                 aa  = hom(k+1,1,30,0)  * ( particles(n)%z / &
    2172                                          ( 0.5 * ( zu(k+1) - zu(k) ) ) )
    2173                 bb  = hom(k+1,1,31,0)  * ( particles(n)%z / &
    2174                                          ( 0.5 * ( zu(k+1) - zu(k) ) ) )
    2175                 cc  = hom(kw+1,1,32,0) * ( particles(n)%z / &
    2176                                          ( 1.0 * ( zw(kw+1) - zw(kw) ) ) )
    2177              ELSE
    2178                 aa  = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
    2179                            ( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
    2180                 bb  = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
    2181                            ( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
    2182                 cc  = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *&
    2183                            ( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
    2184              ENDIF
    2185 
    2186              vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc )
    2187 
    2188              fs_int = ( 2.0 / 3.0 ) * e_mean_int / &
    2189                          ( vv_int + ( 2.0 / 3.0 ) * e_mean_int )
    2190 
    2191 !
    2192 !--          Calculate the Lagrangian timescale according to the suggestion of
    2193 !--          Weil et al. (2004).
    2194              lagr_timescale = ( 4.0 * e_int ) / &
    2195                               ( 3.0 * fs_int * c_0 * diss_int )
    2196 
    2197 !
    2198 !--          Calculate the next particle timestep. dt_gap is the time needed to
    2199 !--          complete the current LES timestep.
    2200              dt_gap = dt_3d - particles(n)%dt_sum
    2201              dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap )
    2202 
    2203 !
    2204 !--          The particle timestep should not be too small in order to prevent
    2205 !--          the number of particle timesteps of getting too large
    2206              IF ( dt_particle < dt_min_part   .AND.  dt_min_part < dt_gap ) &
    2207              THEN
    2208                 dt_particle = dt_min_part
    2209              ENDIF
    2210 
    2211 !
    2212 !--          Calculate the SGS velocity components
    2213              IF ( particles(n)%age == 0.0 )  THEN
    2214 !
    2215 !--             For new particles the SGS components are derived from the SGS
    2216 !--             TKE. Limit the Gaussian random number to the interval
    2217 !--             [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
    2218 !--             from becoming unrealistically large.
    2219                 particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) *   &
    2220                                            ( random_gauss( iran_part, 5.0 ) &
    2221                                              - 1.0 )
    2222                 particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) *   &
    2223                                            ( random_gauss( iran_part, 5.0 ) &
    2224                                              - 1.0 )
    2225                 particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) *   &
    2226                                            ( random_gauss( iran_part, 5.0 ) &
    2227                                              - 1.0 )
    2228 
    2229              ELSE
    2230 
    2231 !
    2232 !--             Restriction of the size of the new timestep: compared to the
    2233 !--             previous timestep the increase must not exceed 200%
    2234 
    2235                 dt_particle_m = particles(n)%age - particles(n)%age_m
    2236                 IF ( dt_particle > 2.0 * dt_particle_m )  THEN
    2237                    dt_particle = 2.0 * dt_particle_m
    2238                 ENDIF
    2239 
    2240 !
    2241 !--             For old particles the SGS components are correlated with the
    2242 !--             values from the previous timestep. Random numbers have also to
    2243 !--             be limited (see above).
    2244 !--             As negative values for the subgrid TKE are not allowed, the
    2245 !--             change of the subgrid TKE with time cannot be smaller than
    2246 !--             -e_int/dt_particle. This value is used as a lower boundary
    2247 !--             value for the change of TKE
    2248 
    2249                 de_dt_min = - e_int / dt_particle
    2250 
    2251                 de_dt = ( e_int - particles(n)%e_m ) / dt_particle_m
    2252 
    2253                 IF ( de_dt < de_dt_min )  THEN
    2254                    de_dt = de_dt_min
    2255                 ENDIF
    2256 
    2257                 particles(n)%rvar1 = particles(n)%rvar1 -                      &
    2258                        fs_int * c_0 * diss_int * particles(n)%rvar1 *          &
    2259                        dt_particle / ( 4.0 * sgs_wfu_part * e_int ) +          &
    2260                        ( 2.0 * sgs_wfu_part * de_dt *                          &
    2261                                particles(n)%rvar1 /                            &
    2262                          ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int            &
    2263                        ) * dt_particle / 2.0                        +          &
    2264                        SQRT( fs_int * c_0 * diss_int ) *                       &
    2265                        ( random_gauss( iran_part, 5.0 ) - 1.0 ) *              &
    2266                        SQRT( dt_particle )
    2267 
    2268                 particles(n)%rvar2 = particles(n)%rvar2 -                      &
    2269                        fs_int * c_0 * diss_int * particles(n)%rvar2 *          &
    2270                        dt_particle / ( 4.0 * sgs_wfv_part * e_int ) +          &
    2271                        ( 2.0 * sgs_wfv_part * de_dt *                          &
    2272                                particles(n)%rvar2 /                            &
    2273                          ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int            &
    2274                        ) * dt_particle / 2.0                        +          &
    2275                        SQRT( fs_int * c_0 * diss_int ) *                       &
    2276                        ( random_gauss( iran_part, 5.0 ) - 1.0 ) *              &
    2277                        SQRT( dt_particle )
    2278 
    2279                 particles(n)%rvar3 = particles(n)%rvar3 -          &
    2280                        fs_int * c_0 * diss_int * particles(n)%rvar3 *    &
    2281                        dt_particle / ( 4.0 * sgs_wfw_part * e_int ) +          &
    2282                        ( 2.0 * sgs_wfw_part * de_dt *                          &
    2283                                particles(n)%rvar3 /                      &
    2284                          ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int            &
    2285                        ) * dt_particle / 2.0                        +          &
    2286                        SQRT( fs_int * c_0 * diss_int ) *                       &
    2287                        ( random_gauss( iran_part, 5.0 ) - 1.0 ) *              &
    2288                        SQRT( dt_particle )
    2289 
    2290              ENDIF
    2291 
    2292              u_int = u_int + particles(n)%rvar1
    2293              v_int = v_int + particles(n)%rvar2
    2294              w_int = w_int + particles(n)%rvar3
    2295 
    2296 !
    2297 !--          Store the SGS TKE of the current timelevel which is needed for
    2298 !--          for calculating the SGS particle velocities at the next timestep
    2299              particles(n)%e_m = e_int
    2300 
    2301           ELSE
    2302 !
    2303 !--          If no SGS velocities are used, only the particle timestep has to
    2304 !--          be set
    2305              dt_particle = dt_3d
    2306 
    2307           ENDIF
    2308 
    2309 !
    2310 !--       Remember the old age of the particle ( needed to prevent that a
    2311 !--       particle crosses several PEs during one timestep and for the
    2312 !--       evaluation of the subgrid particle velocity fluctuations )
    2313           particles(n)%age_m = particles(n)%age
    2314 
    2315 
    2316 !
    2317 !--       Particle advection
    2318           IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 )  THEN
    2319 !
    2320 !--          Pure passive transport (without particle inertia)
    2321              particles(n)%x = particles(n)%x + u_int * dt_particle
    2322              particles(n)%y = particles(n)%y + v_int * dt_particle
    2323              particles(n)%z = particles(n)%z + w_int * dt_particle
    2324 
    2325              particles(n)%speed_x = u_int
    2326              particles(n)%speed_y = v_int
    2327              particles(n)%speed_z = w_int
    2328 
    2329           ELSE
    2330 !
    2331 !--          Transport of particles with inertia
    2332              particles(n)%x = particles(n)%x + particles(n)%speed_x * &
    2333                                                dt_particle
    2334              particles(n)%y = particles(n)%y + particles(n)%speed_y * &
    2335                                                dt_particle
    2336              particles(n)%z = particles(n)%z + particles(n)%speed_z * &
    2337                                                dt_particle
    2338 
    2339 !
    2340 !--          Update of the particle velocity
    2341              dens_ratio = particle_groups(particles(n)%group)%density_ratio
    2342              IF ( cloud_droplets )  THEN
    2343                 exp_arg  = 4.5 * dens_ratio * molecular_viscosity /        &
    2344                            ( particles(n)%radius )**2 *                    &
    2345                            ( 1.0 + 0.15 * ( 2.0 * particles(n)%radius *    &
    2346                              SQRT( ( u_int - particles(n)%speed_x )**2 +   &
    2347                                    ( v_int - particles(n)%speed_y )**2 +   &
    2348                                    ( w_int - particles(n)%speed_z )**2 ) / &
    2349                                             molecular_viscosity )**0.687   &
    2350                            )
    2351                 exp_term = EXP( -exp_arg * dt_particle )
    2352              ELSEIF ( use_sgs_for_particles )  THEN
    2353                 exp_arg  = particle_groups(particles(n)%group)%exp_arg
    2354                 exp_term = EXP( -exp_arg * dt_particle )
    2355              ELSE
    2356                 exp_arg  = particle_groups(particles(n)%group)%exp_arg
    2357                 exp_term = particle_groups(particles(n)%group)%exp_term
    2358              ENDIF
    2359              particles(n)%speed_x = particles(n)%speed_x * exp_term + &
    2360                                     u_int * ( 1.0 - exp_term )
    2361              particles(n)%speed_y = particles(n)%speed_y * exp_term + &
    2362                                     v_int * ( 1.0 - exp_term )
    2363              particles(n)%speed_z = particles(n)%speed_z * exp_term +       &
    2364                              ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) &
    2365                            * ( 1.0 - exp_term )
    2366           ENDIF
    2367 
    2368 !
    2369 !--       Increment the particle age and the total time that the particle
    2370 !--       has advanced within the particle timestep procedure
    2371           particles(n)%age    = particles(n)%age    + dt_particle
    2372           particles(n)%dt_sum = particles(n)%dt_sum + dt_particle
    2373 
    2374 !
    2375 !--       Check whether there is still a particle that has not yet completed
    2376 !--       the total LES timestep
    2377           IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 )  THEN
    2378              dt_3d_reached_l = .FALSE.
    2379           ENDIF
    2380 
    2381        ENDDO   ! advection loop
     178!
     179!--    Particle advection
     180       CALL lpm_advec
    2382181
    2383182!
    2384183!--    Particle reflection from walls
    2385        CALL particle_boundary_conds
     184       CALL lpm_boundary_conds( 'walls' )
    2386185
    2387186!
    2388187!--    User-defined actions after the calculation of the new particle position
    2389        CALL user_advec_particles
     188       CALL user_lpm_advec
    2390189
    2391190!
     
    2400199#endif
    2401200
    2402        CALL cpu_log( log_point_s(44), 'advec_part_advec', 'stop' )
     201       CALL cpu_log( log_point_s(44), 'lpm_advec', 'stop' )
    2403202
    2404203!
    2405204!--    Increment time since last release
    2406205       IF ( dt_3d_reached )  time_prel = time_prel + dt_3d
    2407 
    2408 !    WRITE ( 9, * ) '*** advec_particles: ##0.4'
    2409 !    CALL local_flush( 9 )
    2410 !    nd = 0
    2411 !    DO  n = 1, number_of_particles
    2412 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    2413 !    ENDDO
    2414 !    IF ( nd /= deleted_particles ) THEN
    2415 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    2416 !       CALL local_flush( 9 )
    2417 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    2418 !    ENDIF
    2419206
    2420207!
     
    2423210            dt_3d_reached )  THEN
    2424211
    2425 !
    2426 !--       Check, if particle storage must be extended
    2427           IF ( number_of_particles + number_of_initial_particles > &
    2428                maximum_number_of_particles  )  THEN
    2429              IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    2430                 message_string = 'maximum_number_of_particles ' //   &
    2431                                  'needs to be increased ' //         &
    2432                                  '&but this is not allowed with ' // &
    2433                                  'netcdf_data_format < 3'
    2434                 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
    2435              ELSE
    2436 !    WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory dt_prel'
    2437 !    CALL local_flush( 9 )
    2438                 CALL allocate_prt_memory( number_of_initial_particles )
    2439 !    WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory dt_prel'
    2440 !    CALL local_flush( 9 )
    2441              ENDIF
    2442           ENDIF
    2443 
    2444 !
    2445 !--       Check, if tail storage must be extended
    2446           IF ( use_particle_tails )  THEN
    2447              IF ( number_of_tails + number_of_initial_tails > &
    2448                   maximum_number_of_tails  )  THEN
    2449                 IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    2450                    message_string = 'maximum_number_of_tails ' //    &
    2451                                     'needs to be increased ' //      &
    2452                                     '&but this is not allowed wi' // &
    2453                                     'th netcdf_data_format < 3'
    2454                    CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
    2455                 ELSE
    2456 !    WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory dt_prel'
    2457 !    CALL local_flush( 9 )
    2458                    CALL allocate_tail_memory( number_of_initial_tails )
    2459 !    WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory dt_prel'
    2460 !    CALL local_flush( 9 )
    2461                 ENDIF
    2462              ENDIF
    2463           ENDIF
     212          CALL lpm_release_set
    2464213
    2465214!
     
    2467216!--       restart runs.
    2468217          time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) )
    2469           IF ( number_of_initial_particles /= 0 )  THEN
    2470              is = number_of_particles+1
    2471              ie = number_of_particles+number_of_initial_particles
    2472              particles(is:ie) = initial_particles(1:number_of_initial_particles)
    2473 !
    2474 !--          Add random fluctuation to particle positions. Particles should
    2475 !--          remain in the subdomain.
    2476              IF ( random_start_position )  THEN
    2477                 DO  n = is, ie
    2478                    IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) &
    2479                    THEN
    2480                       particles(n)%x = particles(n)%x + &
    2481                                        ( random_function( iran_part ) - 0.5 ) *&
    2482                                        pdx(particles(n)%group)
    2483                       IF ( particles(n)%x  <=  ( nxl - 0.5 ) * dx )  THEN
    2484                          particles(n)%x = ( nxl - 0.4999999999 ) * dx
    2485                       ELSEIF ( particles(n)%x  >=  ( nxr + 0.5 ) * dx )  THEN
    2486                          particles(n)%x = ( nxr + 0.4999999999 ) * dx
    2487                       ENDIF
    2488                    ENDIF
    2489                    IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) &
    2490                    THEN
    2491                       particles(n)%y = particles(n)%y + &
    2492                                        ( random_function( iran_part ) - 0.5 ) *&
    2493                                        pdy(particles(n)%group)
    2494                       IF ( particles(n)%y  <=  ( nys - 0.5 ) * dy )  THEN
    2495                          particles(n)%y = ( nys - 0.4999999999 ) * dy
    2496                       ELSEIF ( particles(n)%y  >=  ( nyn + 0.5 ) * dy )  THEN
    2497                          particles(n)%y = ( nyn + 0.4999999999 ) * dy
    2498                       ENDIF
    2499                    ENDIF
    2500                    IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) &
    2501                    THEN
    2502                       particles(n)%z = particles(n)%z + &
    2503                                        ( random_function( iran_part ) - 0.5 ) *&
    2504                                        pdz(particles(n)%group)
    2505                    ENDIF
    2506                 ENDDO
    2507              ENDIF
    2508 
    2509 !
    2510 !--          Set the beginning of the new particle tails and their age
    2511              IF ( use_particle_tails )  THEN
    2512                 DO  n = is, ie
    2513 !
    2514 !--                New particles which should have a tail, already have got a
    2515 !--                provisional tail id unequal zero (see init_particles)
    2516                    IF ( particles(n)%tail_id /= 0 )  THEN
    2517                       number_of_tails = number_of_tails + 1
    2518                       nn = number_of_tails
    2519                       particles(n)%tail_id = nn   ! set the final tail id
    2520                       particle_tail_coordinates(1,1,nn) = particles(n)%x
    2521                       particle_tail_coordinates(1,2,nn) = particles(n)%y
    2522                       particle_tail_coordinates(1,3,nn) = particles(n)%z
    2523                       particle_tail_coordinates(1,4,nn) = particles(n)%class
    2524                       particles(n)%tailpoints = 1
    2525                       IF ( minimum_tailpoint_distance /= 0.0 )  THEN
    2526                          particle_tail_coordinates(2,1,nn) = particles(n)%x
    2527                          particle_tail_coordinates(2,2,nn) = particles(n)%y
    2528                          particle_tail_coordinates(2,3,nn) = particles(n)%z
    2529                          particle_tail_coordinates(2,4,nn) = particles(n)%class
    2530                          particle_tail_coordinates(1:2,5,nn) = 0.0
    2531                          particles(n)%tailpoints = 2
    2532                       ENDIF
    2533                    ENDIF
    2534                 ENDDO
    2535              ENDIF
    2536 !    WRITE ( 9, * ) '*** advec_particles: after setting the beginning of new tails'
    2537 !    CALL local_flush( 9 )
    2538 
    2539              number_of_particles = number_of_particles + &
    2540                                    number_of_initial_particles
    2541           ENDIF
    2542218
    2543219       ENDIF
    2544220
    2545 !    WRITE ( 9, * ) '*** advec_particles: ##0.5'
    2546 !    CALL local_flush( 9 )
    2547 !    nd = 0
    2548 !    DO  n = 1, number_of_particles
    2549 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    2550 !    ENDDO
    2551 !    IF ( nd /= deleted_particles ) THEN
    2552 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    2553 !       CALL local_flush( 9 )
    2554 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    2555 !    ENDIF
    2556 
    2557 !    IF ( number_of_particles /= number_of_tails )  THEN
    2558 !       WRITE (9,*) '--- advec_particles: #2'
    2559 !       WRITE (9,*) '    #of p=',number_of_particles,' #of t=',number_of_tails
    2560 !       CALL local_flush( 9 )
    2561 !    ENDIF
    2562 !    DO  n = 1, number_of_particles
    2563 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    2564 !       THEN
    2565 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    2566 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    2567 !          CALL local_flush( 9 )
    2568 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    2569 !       ENDIF
    2570 !    ENDDO
    2571 
    2572 #if defined( __parallel )
    2573 !
    2574 !--    As soon as one particle has moved beyond the boundary of the domain, it
    2575 !--    is included in the relevant transfer arrays and marked for subsequent
    2576 !--    deletion on this PE.
    2577 !--    First run for crossings in x direction. Find out first the number of
    2578 !--    particles to be transferred and allocate temporary arrays needed to store
    2579 !--    them.
    2580 !--    For a one-dimensional decomposition along y, no transfer is necessary,
    2581 !--    because the particle remains on the PE.
    2582        trlp_count  = 0
    2583        trlpt_count = 0
    2584        trrp_count  = 0
    2585        trrpt_count = 0
    2586        IF ( pdims(1) /= 1 )  THEN
    2587 !
    2588 !--       First calculate the storage necessary for sending and receiving the
    2589 !--       data
    2590           DO  n = 1, number_of_particles
    2591              i = ( particles(n)%x + 0.5 * dx ) * ddx
    2592 !
    2593 !--          Above calculation does not work for indices less than zero
    2594              IF ( particles(n)%x < -0.5 * dx )  i = -1
    2595 
    2596              IF ( i < nxl )  THEN
    2597                 trlp_count = trlp_count + 1
    2598                 IF ( particles(n)%tail_id /= 0 )  trlpt_count = trlpt_count + 1
    2599              ELSEIF ( i > nxr )  THEN
    2600                 trrp_count = trrp_count + 1
    2601                 IF ( particles(n)%tail_id /= 0 )  trrpt_count = trrpt_count + 1
    2602              ENDIF
    2603           ENDDO
    2604           IF ( trlp_count  == 0 )  trlp_count  = 1
    2605           IF ( trlpt_count == 0 )  trlpt_count = 1
    2606           IF ( trrp_count  == 0 )  trrp_count  = 1
    2607           IF ( trrpt_count == 0 )  trrpt_count = 1
    2608 
    2609           ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
    2610 
    2611           trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    2612                                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    2613                                 0.0, 0, 0, 0, 0 )
    2614           trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    2615                                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    2616                                 0.0, 0, 0, 0, 0 )
    2617 
    2618           IF ( use_particle_tails )  THEN
    2619              ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
    2620                        trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
    2621              tlength = maximum_number_of_tailpoints * 5
    2622           ENDIF
    2623 
    2624           trlp_count  = 0
    2625           trlpt_count = 0
    2626           trrp_count  = 0
    2627           trrpt_count = 0
    2628 
    2629        ENDIF
    2630 
    2631 !    WRITE ( 9, * ) '*** advec_particles: ##1'
    2632 !    CALL local_flush( 9 )
    2633 !    nd = 0
    2634 !    DO  n = 1, number_of_particles
    2635 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    2636 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    2637 !       THEN
    2638 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    2639 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    2640 !          CALL local_flush( 9 )
    2641 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    2642 !       ENDIF
    2643 !    ENDDO
    2644 !    IF ( nd /= deleted_particles ) THEN
    2645 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    2646 !       CALL local_flush( 9 )
    2647 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    2648 !    ENDIF
    2649 
    2650        DO  n = 1, number_of_particles
    2651 
    2652           nn = particles(n)%tail_id
    2653 
    2654           i = ( particles(n)%x + 0.5 * dx ) * ddx
    2655 !
    2656 !--       Above calculation does not work for indices less than zero
    2657           IF ( particles(n)%x < - 0.5 * dx )  i = -1
    2658 
    2659           IF ( i <  nxl )  THEN
    2660              IF ( i < 0 )  THEN
    2661 !
    2662 !--             Apply boundary condition along x
    2663                 IF ( ibc_par_lr == 0 )  THEN
    2664 !
    2665 !--                Cyclic condition
    2666                    IF ( pdims(1) == 1 )  THEN
    2667                       particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
    2668                       particles(n)%origin_x = ( nx + 1 ) * dx + &
    2669                                               particles(n)%origin_x
    2670                       IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    2671                          i  = particles(n)%tailpoints
    2672                          particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
    2673                                            + particle_tail_coordinates(1:i,1,nn)
    2674                       ENDIF
    2675                    ELSE
    2676                       trlp_count = trlp_count + 1
    2677                       trlp(trlp_count) = particles(n)
    2678                       trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
    2679                       trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
    2680                                                   ( nx + 1 ) * dx
    2681                       particle_mask(n)  = .FALSE.
    2682                       deleted_particles = deleted_particles + 1
    2683 
    2684                       IF ( trlp(trlp_count)%x >= (nx + 0.5)* dx - 1.e-12 ) THEN
    2685                          trlp(trlp_count)%x = trlp(trlp_count)%x - 1.e-10
    2686                          trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x &
    2687                                                      - 1
    2688                       ENDIF
    2689 
    2690                       IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    2691                          trlpt_count = trlpt_count + 1
    2692                          trlpt(:,:,trlpt_count) = &
    2693                                                particle_tail_coordinates(:,:,nn)
    2694                          trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
    2695                                                   trlpt(:,1,trlpt_count)
    2696                          tail_mask(nn) = .FALSE.
    2697                          deleted_tails = deleted_tails + 1
    2698                       ENDIF
    2699                    ENDIF
    2700 
    2701                 ELSEIF ( ibc_par_lr == 1 )  THEN
    2702 !
    2703 !--                Particle absorption
    2704                    particle_mask(n) = .FALSE.
    2705                    deleted_particles = deleted_particles + 1
    2706                    IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    2707                       tail_mask(nn) = .FALSE.
    2708                       deleted_tails = deleted_tails + 1
    2709                    ENDIF
    2710 
    2711                 ELSEIF ( ibc_par_lr == 2 )  THEN
    2712 !
    2713 !--                Particle reflection
    2714                    particles(n)%x       = -particles(n)%x
    2715                    particles(n)%speed_x = -particles(n)%speed_x
    2716 
    2717                 ENDIF
    2718              ELSE
    2719 !
    2720 !--             Store particle data in the transfer array, which will be send
    2721 !--             to the neighbouring PE
    2722                 trlp_count = trlp_count + 1
    2723                 trlp(trlp_count) = particles(n)
    2724                 particle_mask(n) = .FALSE.
    2725                 deleted_particles = deleted_particles + 1
    2726 
    2727                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    2728                    trlpt_count = trlpt_count + 1
    2729                    trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
    2730                    tail_mask(nn) = .FALSE.
    2731                    deleted_tails = deleted_tails + 1
    2732                 ENDIF
    2733             ENDIF
    2734 
    2735           ELSEIF ( i > nxr )  THEN
    2736              IF ( i > nx )  THEN
    2737 !
    2738 !--             Apply boundary condition along x
    2739                 IF ( ibc_par_lr == 0 )  THEN
    2740 !
    2741 !--                Cyclic condition
    2742                    IF ( pdims(1) == 1 )  THEN
    2743                       particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
    2744                       particles(n)%origin_x = particles(n)%origin_x - &
    2745                                               ( nx + 1 ) * dx
    2746                       IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    2747                          i = particles(n)%tailpoints
    2748                          particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
    2749                                            + particle_tail_coordinates(1:i,1,nn)
    2750                       ENDIF
    2751                    ELSE
    2752                       trrp_count = trrp_count + 1
    2753                       trrp(trrp_count) = particles(n)
    2754                       trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
    2755                       trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
    2756                                                   ( nx + 1 ) * dx
    2757                       particle_mask(n) = .FALSE.
    2758                       deleted_particles = deleted_particles + 1
    2759 
    2760                       IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    2761                          trrpt_count = trrpt_count + 1
    2762                          trrpt(:,:,trrpt_count) = &
    2763                                                particle_tail_coordinates(:,:,nn)
    2764                          trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
    2765                                                   ( nx + 1 ) * dx
    2766                          tail_mask(nn) = .FALSE.
    2767                          deleted_tails = deleted_tails + 1
    2768                       ENDIF
    2769                    ENDIF
    2770 
    2771                 ELSEIF ( ibc_par_lr == 1 )  THEN
    2772 !
    2773 !--                Particle absorption
    2774                    particle_mask(n) = .FALSE.
    2775                    deleted_particles = deleted_particles + 1
    2776                    IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    2777                       tail_mask(nn) = .FALSE.
    2778                       deleted_tails = deleted_tails + 1
    2779                    ENDIF
    2780 
    2781                 ELSEIF ( ibc_par_lr == 2 )  THEN
    2782 !
    2783 !--                Particle reflection
    2784                    particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
    2785                    particles(n)%speed_x = -particles(n)%speed_x
    2786 
    2787                 ENDIF
    2788              ELSE
    2789 !
    2790 !--             Store particle data in the transfer array, which will be send
    2791 !--             to the neighbouring PE
    2792                 trrp_count = trrp_count + 1
    2793                 trrp(trrp_count) = particles(n)
    2794                 particle_mask(n) = .FALSE.
    2795                 deleted_particles = deleted_particles + 1
    2796 
    2797                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    2798                    trrpt_count = trrpt_count + 1
    2799                    trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
    2800                    tail_mask(nn) = .FALSE.
    2801                    deleted_tails = deleted_tails + 1
    2802                 ENDIF
    2803              ENDIF
    2804 
    2805           ENDIF
    2806        ENDDO
    2807 
    2808 !    WRITE ( 9, * ) '*** advec_particles: ##2'
    2809 !    CALL local_flush( 9 )
    2810 !    nd = 0
    2811 !    DO  n = 1, number_of_particles
    2812 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    2813 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    2814 !       THEN
    2815 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    2816 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    2817 !          CALL local_flush( 9 )
    2818 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    2819 !       ENDIF
    2820 !    ENDDO
    2821 !    IF ( nd /= deleted_particles ) THEN
    2822 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    2823 !       CALL local_flush( 9 )
    2824 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    2825 !    ENDIF
    2826 
    2827 !
    2828 !--    Send left boundary, receive right boundary (but first exchange how many
    2829 !--    and check, if particle storage must be extended)
    2830        IF ( pdims(1) /= 1 )  THEN
    2831 
    2832           CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'start' )
    2833           CALL MPI_SENDRECV( trlp_count,      1, MPI_INTEGER, pleft,  0, &
    2834                              trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
    2835                              comm2d, status, ierr )
    2836 
    2837           IF ( number_of_particles + trrp_count_recv > &
    2838                maximum_number_of_particles )           &
    2839           THEN
    2840              IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    2841                  message_string = 'maximum_number_of_particles ' //    &
    2842                                   'needs to be increased ' //          &
    2843                                   '&but this is not allowed with ' //  &
    2844                                   'netcdf-data_format < 3'
    2845                 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
    2846              ELSE
    2847 !    WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trrp'
    2848 !    CALL local_flush( 9 )
    2849                 CALL allocate_prt_memory( trrp_count_recv )
    2850 !    WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trrp'
    2851 !    CALL local_flush( 9 )
    2852              ENDIF
    2853           ENDIF
    2854 
    2855           CALL MPI_SENDRECV( trlp(1)%age, trlp_count, mpi_particle_type,     &
    2856                              pleft, 1, particles(number_of_particles+1)%age, &
    2857                              trrp_count_recv, mpi_particle_type, pright, 1,  &
    2858                              comm2d, status, ierr )
    2859 
    2860           IF ( use_particle_tails )  THEN
    2861 
    2862              CALL MPI_SENDRECV( trlpt_count,      1, MPI_INTEGER, pleft,  0, &
    2863                                 trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
    2864                                 comm2d, status, ierr )
    2865 
    2866              IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
    2867              THEN
    2868                 IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    2869                    message_string = 'maximum_number_of_tails ' //   &
    2870                                     'needs to be increased ' //     &
    2871                                     '&but this is not allowed wi'// &
    2872                                     'th netcdf_data_format < 3'
    2873                    CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
    2874                 ELSE
    2875 !    WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trrpt'
    2876 !    CALL local_flush( 9 )
    2877                    CALL allocate_tail_memory( trrpt_count_recv )
    2878 !    WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trrpt'
    2879 !    CALL local_flush( 9 )
    2880                 ENDIF
    2881              ENDIF
    2882 
    2883              CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL,   &
    2884                                 pleft, 1,                                      &
    2885                              particle_tail_coordinates(1,1,number_of_tails+1), &
    2886                                 trrpt_count_recv*tlength, MPI_REAL, pright, 1, &
    2887                                 comm2d, status, ierr )
    2888 !
    2889 !--          Update the tail ids for the transferred particles
    2890              nn = number_of_tails
    2891              DO  n = number_of_particles+1, number_of_particles+trrp_count_recv
    2892                 IF ( particles(n)%tail_id /= 0 )  THEN
    2893                    nn = nn + 1
    2894                    particles(n)%tail_id = nn
    2895                 ENDIF
    2896              ENDDO
    2897 
    2898           ENDIF
    2899 
    2900           number_of_particles = number_of_particles + trrp_count_recv
    2901           number_of_tails     = number_of_tails     + trrpt_count_recv
    2902 !       IF ( number_of_particles /= number_of_tails )  THEN
    2903 !          WRITE (9,*) '--- advec_particles: #3'
    2904 !          WRITE (9,*) '    #of p=',number_of_particles,' #of t=',number_of_tails
    2905 !          CALL local_flush( 9 )
    2906 !       ENDIF
    2907 
    2908 !
    2909 !--       Send right boundary, receive left boundary
    2910           CALL MPI_SENDRECV( trrp_count,      1, MPI_INTEGER, pright, 0, &
    2911                              trlp_count_recv, 1, MPI_INTEGER, pleft,  0, &
    2912                              comm2d, status, ierr )
    2913 
    2914           IF ( number_of_particles + trlp_count_recv > &
    2915                maximum_number_of_particles )           &
    2916           THEN
    2917              IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    2918                 message_string = 'maximum_number_of_particles ' //  &
    2919                                  'needs to be increased ' //        &
    2920                                  '&but this is not allowed with '// &
    2921                                  'netcdf_data_format < 3'
    2922                 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
    2923              ELSE
    2924 !    WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trlp'
    2925 !    CALL local_flush( 9 )
    2926                 CALL allocate_prt_memory( trlp_count_recv )
    2927 !    WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trlp'
    2928 !    CALL local_flush( 9 )
    2929              ENDIF
    2930           ENDIF
    2931 
    2932           CALL MPI_SENDRECV( trrp(1)%age, trrp_count, mpi_particle_type,      &
    2933                              pright, 1, particles(number_of_particles+1)%age, &
    2934                              trlp_count_recv, mpi_particle_type, pleft, 1,    &
    2935                              comm2d, status, ierr )
    2936 
    2937           IF ( use_particle_tails )  THEN
    2938 
    2939              CALL MPI_SENDRECV( trrpt_count,      1, MPI_INTEGER, pright, 0, &
    2940                                 trlpt_count_recv, 1, MPI_INTEGER, pleft,  0, &
    2941                                 comm2d, status, ierr )
    2942 
    2943              IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
    2944              THEN
    2945                 IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    2946                    message_string = 'maximum_number_of_tails ' //   &
    2947                                     'needs to be increased ' //     &
    2948                                     '&but this is not allowed wi'// &
    2949                                     'th netcdf_data_format < 3'
    2950                    CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
    2951                 ELSE
    2952 !    WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trlpt'
    2953 !    CALL local_flush( 9 )
    2954                    CALL allocate_tail_memory( trlpt_count_recv )
    2955 !    WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trlpt'
    2956 !    CALL local_flush( 9 )
    2957                 ENDIF
    2958              ENDIF
    2959 
    2960              CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL,   &
    2961                                 pright, 1,                                     &
    2962                              particle_tail_coordinates(1,1,number_of_tails+1), &
    2963                                 trlpt_count_recv*tlength, MPI_REAL, pleft, 1,  &
    2964                                 comm2d, status, ierr )
    2965 !
    2966 !--          Update the tail ids for the transferred particles
    2967              nn = number_of_tails
    2968              DO  n = number_of_particles+1, number_of_particles+trlp_count_recv
    2969                 IF ( particles(n)%tail_id /= 0 )  THEN
    2970                    nn = nn + 1
    2971                    particles(n)%tail_id = nn
    2972                 ENDIF
    2973              ENDDO
    2974 
    2975           ENDIF
    2976 
    2977           number_of_particles = number_of_particles + trlp_count_recv
    2978           number_of_tails     = number_of_tails     + trlpt_count_recv
    2979 !       IF ( number_of_particles /= number_of_tails )  THEN
    2980 !          WRITE (9,*) '--- advec_particles: #4'
    2981 !          WRITE (9,*) '    #of p=',number_of_particles,' #of t=',number_of_tails
    2982 !          CALL local_flush( 9 )
    2983 !       ENDIF
    2984 
    2985           IF ( use_particle_tails )  THEN
    2986              DEALLOCATE( trlpt, trrpt )
    2987           ENDIF
    2988           DEALLOCATE( trlp, trrp )
    2989 
    2990           CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'pause' )
    2991 
    2992        ENDIF
    2993 
    2994 !    WRITE ( 9, * ) '*** advec_particles: ##3'
    2995 !    CALL local_flush( 9 )
    2996 !    nd = 0
    2997 !    DO  n = 1, number_of_particles
    2998 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    2999 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    3000 !       THEN
    3001 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3002 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    3003 !          CALL local_flush( 9 )
    3004 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3005 !       ENDIF
    3006 !    ENDDO
    3007 !    IF ( nd /= deleted_particles ) THEN
    3008 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    3009 !       CALL local_flush( 9 )
    3010 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    3011 !    ENDIF
    3012 
    3013 !
    3014 !--    Check whether particles have crossed the boundaries in y direction. Note
    3015 !--    that this case can also apply to particles that have just been received
    3016 !--    from the adjacent right or left PE.
    3017 !--    Find out first the number of particles to be transferred and allocate
    3018 !--    temporary arrays needed to store them.
    3019 !--    For a one-dimensional decomposition along x, no transfer is necessary,
    3020 !--    because the particle remains on the PE.
    3021        trsp_count  = 0
    3022        trspt_count = 0
    3023        trnp_count  = 0
    3024        trnpt_count = 0
    3025        IF ( pdims(2) /= 1 )  THEN
    3026 !
    3027 !--       First calculate the storage necessary for sending and receiving the
    3028 !--       data
    3029           DO  n = 1, number_of_particles
    3030              IF ( particle_mask(n) )  THEN
    3031                 j = ( particles(n)%y + 0.5 * dy ) * ddy
    3032 !
    3033 !--             Above calculation does not work for indices less than zero
    3034                 IF ( particles(n)%y < -0.5 * dy )  j = -1
    3035 
    3036                 IF ( j < nys )  THEN
    3037                    trsp_count = trsp_count + 1
    3038                    IF ( particles(n)%tail_id /= 0 )  trspt_count = trspt_count+1
    3039                 ELSEIF ( j > nyn )  THEN
    3040                    trnp_count = trnp_count + 1
    3041                    IF ( particles(n)%tail_id /= 0 )  trnpt_count = trnpt_count+1
    3042                 ENDIF
    3043              ENDIF
    3044           ENDDO
    3045           IF ( trsp_count  == 0 )  trsp_count  = 1
    3046           IF ( trspt_count == 0 )  trspt_count = 1
    3047           IF ( trnp_count  == 0 )  trnp_count  = 1
    3048           IF ( trnpt_count == 0 )  trnpt_count = 1
    3049 
    3050           ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
    3051 
    3052           trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    3053                                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    3054                                 0.0, 0, 0, 0, 0 )
    3055           trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    3056                                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    3057                                 0.0, 0, 0, 0, 0 )
    3058 
    3059           IF ( use_particle_tails )  THEN
    3060              ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
    3061                        trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
    3062              tlength = maximum_number_of_tailpoints * 5
    3063           ENDIF
    3064 
    3065           trsp_count  = 0
    3066           trspt_count = 0
    3067           trnp_count  = 0
    3068           trnpt_count = 0
    3069 
    3070        ENDIF
    3071 
    3072 !    WRITE ( 9, * ) '*** advec_particles: ##4'
    3073 !    CALL local_flush( 9 )
    3074 !    nd = 0
    3075 !    DO  n = 1, number_of_particles
    3076 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    3077 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    3078 !       THEN
    3079 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3080 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    3081 !          CALL local_flush( 9 )
    3082 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3083 !       ENDIF
    3084 !    ENDDO
    3085 !    IF ( nd /= deleted_particles ) THEN
    3086 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    3087 !       CALL local_flush( 9 )
    3088 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    3089 !    ENDIF
    3090 
    3091        DO  n = 1, number_of_particles
    3092 
    3093           nn = particles(n)%tail_id
    3094 !
    3095 !--       Only those particles that have not been marked as 'deleted' may be
    3096 !--       moved.
    3097           IF ( particle_mask(n) )  THEN
    3098              j = ( particles(n)%y + 0.5 * dy ) * ddy
    3099 !
    3100 !--          Above calculation does not work for indices less than zero
    3101              IF ( particles(n)%y < -0.5 * dy )  j = -1
    3102 
    3103              IF ( j < nys )  THEN
    3104                 IF ( j < 0 )  THEN
    3105 !
    3106 !--                Apply boundary condition along y
    3107                    IF ( ibc_par_ns == 0 )  THEN
    3108 !
    3109 !--                   Cyclic condition
    3110                       IF ( pdims(2) == 1 )  THEN
    3111                          particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
    3112                          particles(n)%origin_y = ( ny + 1 ) * dy + &
    3113                                                  particles(n)%origin_y
    3114                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3115                             i = particles(n)%tailpoints
    3116                             particle_tail_coordinates(1:i,2,nn) = ( ny+1 ) * dy&
    3117                                            + particle_tail_coordinates(1:i,2,nn)
    3118                          ENDIF
    3119                       ELSE
    3120                          trsp_count = trsp_count + 1
    3121                          trsp(trsp_count) = particles(n)
    3122                          trsp(trsp_count)%y = ( ny + 1 ) * dy + &
    3123                                               trsp(trsp_count)%y
    3124                          trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
    3125                                               + ( ny + 1 ) * dy
    3126                          particle_mask(n) = .FALSE.
    3127                          deleted_particles = deleted_particles + 1
    3128 
    3129                          IF ( trsp(trsp_count)%y >= (ny+0.5)* dy - 1.e-12 ) THEN
    3130                             trsp(trsp_count)%y = trsp(trsp_count)%y - 1.e-10
    3131                             trsp(trsp_count)%origin_y =                        &
    3132                                                     trsp(trsp_count)%origin_y - 1
    3133                          ENDIF
    3134 
    3135                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3136                             trspt_count = trspt_count + 1
    3137                             trspt(:,:,trspt_count) = &
    3138                                                particle_tail_coordinates(:,:,nn)
    3139                             trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
    3140                                                      trspt(:,2,trspt_count)
    3141                             tail_mask(nn) = .FALSE.
    3142                             deleted_tails = deleted_tails + 1
    3143                          ENDIF
    3144                       ENDIF
    3145 
    3146                    ELSEIF ( ibc_par_ns == 1 )  THEN
    3147 !
    3148 !--                   Particle absorption
    3149                       particle_mask(n) = .FALSE.
    3150                       deleted_particles = deleted_particles + 1
    3151                       IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3152                          tail_mask(nn) = .FALSE.
    3153                          deleted_tails = deleted_tails + 1
    3154                       ENDIF
    3155 
    3156                    ELSEIF ( ibc_par_ns == 2 )  THEN
    3157 !
    3158 !--                   Particle reflection
    3159                       particles(n)%y       = -particles(n)%y
    3160                       particles(n)%speed_y = -particles(n)%speed_y
    3161 
    3162                    ENDIF
    3163                 ELSE
    3164 !
    3165 !--                Store particle data in the transfer array, which will be send
    3166 !--                to the neighbouring PE
    3167                    trsp_count = trsp_count + 1
    3168                    trsp(trsp_count) = particles(n)
    3169                    particle_mask(n) = .FALSE.
    3170                    deleted_particles = deleted_particles + 1
    3171 
    3172                    IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3173                       trspt_count = trspt_count + 1
    3174                       trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
    3175                       tail_mask(nn) = .FALSE.
    3176                       deleted_tails = deleted_tails + 1
    3177                    ENDIF
    3178                 ENDIF
    3179 
    3180              ELSEIF ( j > nyn )  THEN
    3181                 IF ( j > ny )  THEN
    3182 !
    3183 !--                Apply boundary condition along x
    3184                    IF ( ibc_par_ns == 0 )  THEN
    3185 !
    3186 !--                   Cyclic condition
    3187                       IF ( pdims(2) == 1 )  THEN
    3188                          particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
    3189                          particles(n)%origin_y = particles(n)%origin_y - &
    3190                                                  ( ny + 1 ) * dy
    3191                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3192                             i = particles(n)%tailpoints
    3193                             particle_tail_coordinates(1:i,2,nn) = - (ny+1) * dy&
    3194                                            + particle_tail_coordinates(1:i,2,nn)
    3195                          ENDIF
    3196                       ELSE
    3197                          trnp_count = trnp_count + 1
    3198                          trnp(trnp_count) = particles(n)
    3199                          trnp(trnp_count)%y = trnp(trnp_count)%y - &
    3200                                               ( ny + 1 ) * dy
    3201                          trnp(trnp_count)%origin_y = trnp(trnp_count)%origin_y &
    3202                                                      - ( ny + 1 ) * dy
    3203                          particle_mask(n) = .FALSE.
    3204                          deleted_particles = deleted_particles + 1
    3205 
    3206                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3207                             trnpt_count = trnpt_count + 1
    3208                             trnpt(:,:,trnpt_count) = &
    3209                                                particle_tail_coordinates(:,:,nn)
    3210                             trnpt(:,2,trnpt_count) = trnpt(:,2,trnpt_count) - &
    3211                                                      ( ny + 1 ) * dy
    3212                             tail_mask(nn) = .FALSE.
    3213                             deleted_tails = deleted_tails + 1
    3214                          ENDIF
    3215                       ENDIF
    3216 
    3217                    ELSEIF ( ibc_par_ns == 1 )  THEN
    3218 !
    3219 !--                   Particle absorption
    3220                       particle_mask(n) = .FALSE.
    3221                       deleted_particles = deleted_particles + 1
    3222                       IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3223                          tail_mask(nn) = .FALSE.
    3224                          deleted_tails = deleted_tails + 1
    3225                       ENDIF
    3226 
    3227                    ELSEIF ( ibc_par_ns == 2 )  THEN
    3228 !
    3229 !--                   Particle reflection
    3230                       particles(n)%y       = 2 * ( ny * dy ) - particles(n)%y
    3231                       particles(n)%speed_y = -particles(n)%speed_y
    3232 
    3233                    ENDIF
    3234                 ELSE
    3235 !
    3236 !--                Store particle data in the transfer array, which will be send
    3237 !--                to the neighbouring PE
    3238                    trnp_count = trnp_count + 1
    3239                    trnp(trnp_count) = particles(n)
    3240                    particle_mask(n) = .FALSE.
    3241                    deleted_particles = deleted_particles + 1
    3242 
    3243                    IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3244                       trnpt_count = trnpt_count + 1
    3245                       trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
    3246                       tail_mask(nn) = .FALSE.
    3247                       deleted_tails = deleted_tails + 1
    3248                    ENDIF
    3249                 ENDIF
    3250 
    3251              ENDIF
    3252           ENDIF
    3253        ENDDO
    3254 
    3255 !    WRITE ( 9, * ) '*** advec_particles: ##5'
    3256 !    CALL local_flush( 9 )
    3257 !    nd = 0
    3258 !    DO  n = 1, number_of_particles
    3259 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    3260 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    3261 !       THEN
    3262 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3263 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    3264 !          CALL local_flush( 9 )
    3265 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3266 !       ENDIF
    3267 !    ENDDO
    3268 !    IF ( nd /= deleted_particles ) THEN
    3269 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    3270 !       CALL local_flush( 9 )
    3271 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    3272 !    ENDIF
    3273 
    3274 !
    3275 !--    Send front boundary, receive back boundary (but first exchange how many
    3276 !--    and check, if particle storage must be extended)
    3277        IF ( pdims(2) /= 1 )  THEN
    3278 
    3279           CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'continue' )
    3280           CALL MPI_SENDRECV( trsp_count,      1, MPI_INTEGER, psouth, 0, &
    3281                              trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
    3282                              comm2d, status, ierr )
    3283 
    3284           IF ( number_of_particles + trnp_count_recv > &
    3285                maximum_number_of_particles )           &
    3286           THEN
    3287              IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    3288                 message_string = 'maximum_number_of_particles ' //  &
    3289                                  'needs to be increased ' //        &
    3290                                  '&but this is not allowed with '// &
    3291                                  'netcdf_data_format < 3'
    3292                 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
    3293              ELSE
    3294 !    WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trnp'
    3295 !    CALL local_flush( 9 )
    3296                 CALL allocate_prt_memory( trnp_count_recv )
    3297 !    WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trnp'
    3298 !    CALL local_flush( 9 )
    3299              ENDIF
    3300           ENDIF
    3301 
    3302           CALL MPI_SENDRECV( trsp(1)%age, trsp_count, mpi_particle_type,      &
    3303                              psouth, 1, particles(number_of_particles+1)%age, &
    3304                              trnp_count_recv, mpi_particle_type, pnorth, 1,   &
    3305                              comm2d, status, ierr )
    3306 
    3307           IF ( use_particle_tails )  THEN
    3308 
    3309              CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, psouth, 0, &
    3310                                 trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
    3311                                 comm2d, status, ierr )
    3312 
    3313              IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
    3314              THEN
    3315                 IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    3316                    message_string = 'maximum_number_of_tails ' //    &
    3317                                     'needs to be increased ' //      &
    3318                                     '&but this is not allowed wi' // &
    3319                                     'th netcdf_data_format < 3'
    3320                    CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
    3321                 ELSE
    3322 !    WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trnpt'
    3323 !    CALL local_flush( 9 )
    3324                    CALL allocate_tail_memory( trnpt_count_recv )
    3325 !    WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trnpt'
    3326 !    CALL local_flush( 9 )
    3327                 ENDIF
    3328              ENDIF
    3329 
    3330              CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL,   &
    3331                                 psouth, 1,                                     &
    3332                              particle_tail_coordinates(1,1,number_of_tails+1), &
    3333                                 trnpt_count_recv*tlength, MPI_REAL, pnorth, 1, &
    3334                                 comm2d, status, ierr )
    3335 
    3336 !
    3337 !--          Update the tail ids for the transferred particles
    3338              nn = number_of_tails
    3339              DO  n = number_of_particles+1, number_of_particles+trnp_count_recv
    3340                 IF ( particles(n)%tail_id /= 0 )  THEN
    3341                    nn = nn + 1
    3342                    particles(n)%tail_id = nn
    3343                 ENDIF
    3344              ENDDO
    3345 
    3346           ENDIF
    3347 
    3348           number_of_particles = number_of_particles + trnp_count_recv
    3349           number_of_tails     = number_of_tails     + trnpt_count_recv
    3350 !       IF ( number_of_particles /= number_of_tails )  THEN
    3351 !          WRITE (9,*) '--- advec_particles: #5'
    3352 !          WRITE (9,*) '    #of p=',number_of_particles,' #of t=',number_of_tails
    3353 !          CALL local_flush( 9 )
    3354 !       ENDIF
    3355 
    3356 !
    3357 !--       Send back boundary, receive front boundary
    3358           CALL MPI_SENDRECV( trnp_count,      1, MPI_INTEGER, pnorth, 0, &
    3359                              trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
    3360                              comm2d, status, ierr )
    3361 
    3362           IF ( number_of_particles + trsp_count_recv > &
    3363                maximum_number_of_particles )           &
    3364           THEN
    3365              IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    3366                 message_string = 'maximum_number_of_particles ' //   &
    3367                                  'needs to be increased ' //         &
    3368                                  '&but this is not allowed with ' // &
    3369                                  'netcdf_data_format < 3'
    3370                CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 ) 
    3371              ELSE
    3372 !    WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trsp'
    3373 !    CALL local_flush( 9 )
    3374                 CALL allocate_prt_memory( trsp_count_recv )
    3375 !    WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trsp'
    3376 !    CALL local_flush( 9 )
    3377              ENDIF
    3378           ENDIF
    3379 
    3380           CALL MPI_SENDRECV( trnp(1)%age, trnp_count, mpi_particle_type,      &
    3381                              pnorth, 1, particles(number_of_particles+1)%age, &
    3382                              trsp_count_recv, mpi_particle_type, psouth, 1,   &
    3383                              comm2d, status, ierr )
    3384 
    3385           IF ( use_particle_tails )  THEN
    3386 
    3387              CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pnorth, 0, &
    3388                                 trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
    3389                                 comm2d, status, ierr )
    3390 
    3391              IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
    3392              THEN
    3393                 IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
    3394                    message_string = 'maximum_number_of_tails ' //   &
    3395                                     'needs to be increased ' //     &
    3396                                     '&but this is not allowed wi'// &
    3397                                     'th NetCDF output switched on'
    3398                    CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
    3399                 ELSE
    3400 !    WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trspt'
    3401 !    CALL local_flush( 9 )
    3402                    CALL allocate_tail_memory( trspt_count_recv )
    3403 !    WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trspt'
    3404 !    CALL local_flush( 9 )
    3405                 ENDIF
    3406              ENDIF
    3407 
    3408              CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL,   &
    3409                                 pnorth, 1,                                     &
    3410                              particle_tail_coordinates(1,1,number_of_tails+1), &
    3411                                 trspt_count_recv*tlength, MPI_REAL, psouth, 1, &
    3412                                 comm2d, status, ierr )
    3413 !
    3414 !--          Update the tail ids for the transferred particles
    3415              nn = number_of_tails
    3416              DO  n = number_of_particles+1, number_of_particles+trsp_count_recv
    3417                 IF ( particles(n)%tail_id /= 0 )  THEN
    3418                    nn = nn + 1
    3419                    particles(n)%tail_id = nn
    3420                 ENDIF
    3421              ENDDO
    3422 
    3423           ENDIF
    3424 
    3425           number_of_particles = number_of_particles + trsp_count_recv
    3426           number_of_tails     = number_of_tails     + trspt_count_recv
    3427 !       IF ( number_of_particles /= number_of_tails )  THEN
    3428 !          WRITE (9,*) '--- advec_particles: #6'
    3429 !          WRITE (9,*) '    #of p=',number_of_particles,' #of t=',number_of_tails
    3430 !          CALL local_flush( 9 )
    3431 !       ENDIF
    3432 
    3433           IF ( use_particle_tails )  THEN
    3434              DEALLOCATE( trspt, trnpt )
    3435           ENDIF
    3436           DEALLOCATE( trsp, trnp )
    3437 
    3438           CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'stop' )
    3439 
    3440        ENDIF
    3441 
    3442 !    WRITE ( 9, * ) '*** advec_particles: ##6'
    3443 !    CALL local_flush( 9 )
    3444 !    nd = 0
    3445 !    DO  n = 1, number_of_particles
    3446 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    3447 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    3448 !       THEN
    3449 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3450 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    3451 !          CALL local_flush( 9 )
    3452 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3453 !       ENDIF
    3454 !    ENDDO
    3455 !    IF ( nd /= deleted_particles ) THEN
    3456 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    3457 !       CALL local_flush( 9 )
    3458 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    3459 !    ENDIF
    3460 
    3461 #else
    3462 
    3463 !
    3464 !--    Apply boundary conditions
    3465        DO  n = 1, number_of_particles
    3466 
    3467           nn = particles(n)%tail_id
    3468 
    3469           IF ( particles(n)%x < -0.5 * dx )  THEN
    3470 
    3471              IF ( ibc_par_lr == 0 )  THEN
    3472 !
    3473 !--             Cyclic boundary. Relevant coordinate has to be changed.
    3474                 particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
    3475                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3476                    i = particles(n)%tailpoints
    3477                    particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
    3478                                              particle_tail_coordinates(1:i,1,nn)
    3479                 ENDIF
    3480              ELSEIF ( ibc_par_lr == 1 )  THEN
    3481 !
    3482 !--             Particle absorption
    3483                 particle_mask(n) = .FALSE.
    3484                 deleted_particles = deleted_particles + 1
    3485                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3486                    tail_mask(nn) = .FALSE.
    3487                    deleted_tails = deleted_tails + 1
    3488                 ENDIF
    3489              ELSEIF ( ibc_par_lr == 2 )  THEN
    3490 !
    3491 !--             Particle reflection
    3492                 particles(n)%x       = -dx - particles(n)%x
    3493                 particles(n)%speed_x = -particles(n)%speed_x
    3494              ENDIF
    3495 
    3496           ELSEIF ( particles(n)%x >= ( nx + 0.5 ) * dx )  THEN
    3497 
    3498              IF ( ibc_par_lr == 0 )  THEN
    3499 !
    3500 !--             Cyclic boundary. Relevant coordinate has to be changed.
    3501                 particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
    3502                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3503                    i = particles(n)%tailpoints
    3504                    particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
    3505                                              particle_tail_coordinates(1:i,1,nn)
    3506                 ENDIF
    3507              ELSEIF ( ibc_par_lr == 1 )  THEN
    3508 !
    3509 !--             Particle absorption
    3510                 particle_mask(n) = .FALSE.
    3511                 deleted_particles = deleted_particles + 1
    3512                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3513                    tail_mask(nn) = .FALSE.
    3514                    deleted_tails = deleted_tails + 1
    3515                 ENDIF
    3516              ELSEIF ( ibc_par_lr == 2 )  THEN
    3517 !
    3518 !--             Particle reflection
    3519                 particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
    3520                 particles(n)%speed_x = -particles(n)%speed_x
    3521              ENDIF
    3522 
    3523           ENDIF
    3524 
    3525           IF ( particles(n)%y < -0.5 * dy )  THEN
    3526 
    3527              IF ( ibc_par_ns == 0 )  THEN
    3528 !
    3529 !--             Cyclic boundary. Relevant coordinate has to be changed.
    3530                 particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
    3531                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3532                    i = particles(n)%tailpoints
    3533                    particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
    3534                                              particle_tail_coordinates(1:i,2,nn)
    3535                 ENDIF
    3536              ELSEIF ( ibc_par_ns == 1 )  THEN
    3537 !
    3538 !--             Particle absorption
    3539                 particle_mask(n) = .FALSE.
    3540                 deleted_particles = deleted_particles + 1
    3541                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3542                    tail_mask(nn) = .FALSE.
    3543                    deleted_tails = deleted_tails + 1
    3544                 ENDIF
    3545              ELSEIF ( ibc_par_ns == 2 )  THEN
    3546 !
    3547 !--             Particle reflection
    3548                 particles(n)%y       = -dy - particles(n)%y
    3549                 particles(n)%speed_y = -particles(n)%speed_y
    3550              ENDIF
    3551 
    3552           ELSEIF ( particles(n)%y >= ( ny + 0.5 ) * dy )  THEN
    3553 
    3554              IF ( ibc_par_ns == 0 )  THEN
    3555 !
    3556 !--             Cyclic boundary. Relevant coordinate has to be changed.
    3557                 particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
    3558                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3559                    i = particles(n)%tailpoints
    3560                    particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
    3561                                              particle_tail_coordinates(1:i,2,nn)
    3562                 ENDIF
    3563              ELSEIF ( ibc_par_ns == 1 )  THEN
    3564 !
    3565 !--             Particle absorption
    3566                 particle_mask(n) = .FALSE.
    3567                 deleted_particles = deleted_particles + 1
    3568                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3569                    tail_mask(nn) = .FALSE.
    3570                    deleted_tails = deleted_tails + 1
    3571                 ENDIF
    3572              ELSEIF ( ibc_par_ns == 2 )  THEN
    3573 !
    3574 !--             Particle reflection
    3575                 particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
    3576                 particles(n)%speed_y = -particles(n)%speed_y
    3577              ENDIF
    3578 
    3579           ENDIF
    3580        ENDDO
    3581 
    3582 #endif
     221!
     222!--    Horizontal boundary conditions including exchange between subdmains
     223       CALL lpm_exchange_horiz
    3583224
    3584225!
    3585226!--    Apply boundary conditions to those particles that have crossed the top or
    3586227!--    bottom boundary and delete those particles, which are older than allowed
    3587        DO  n = 1, number_of_particles
    3588 
    3589           nn = particles(n)%tail_id
    3590 
    3591 !
    3592 !--       Stop if particles have moved further than the length of one
    3593 !--       PE subdomain (newly released particles have age = age_m!)
    3594           IF ( particles(n)%age /= particles(n)%age_m )  THEN
    3595              IF ( ABS(particles(n)%speed_x) >                                  &
    3596                   ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m)  .OR. &
    3597                   ABS(particles(n)%speed_y) >                                  &
    3598                   ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) )  THEN
    3599 
    3600                   WRITE( message_string, * )  'particle too fast.  n = ',  n
    3601                   CALL message( 'advec_particles', 'PA0148', 2, 2, -1, 6, 1 )
    3602              ENDIF
    3603           ENDIF
    3604 
    3605           IF ( particles(n)%age > particle_maximum_age  .AND.  &
    3606                particle_mask(n) )                              &
    3607           THEN
    3608              particle_mask(n)  = .FALSE.
    3609              deleted_particles = deleted_particles + 1
    3610              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3611                 tail_mask(nn) = .FALSE.
    3612                 deleted_tails = deleted_tails + 1
    3613              ENDIF
    3614           ENDIF
    3615 
    3616           IF ( particles(n)%z >= zu(nz)  .AND.  particle_mask(n) )  THEN
    3617              IF ( ibc_par_t == 1 )  THEN
    3618 !
    3619 !--             Particle absorption
    3620                 particle_mask(n)  = .FALSE.
    3621                 deleted_particles = deleted_particles + 1
    3622                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3623                    tail_mask(nn) = .FALSE.
    3624                    deleted_tails = deleted_tails + 1
    3625                 ENDIF
    3626              ELSEIF ( ibc_par_t == 2 )  THEN
    3627 !
    3628 !--             Particle reflection
    3629                 particles(n)%z       = 2.0 * zu(nz) - particles(n)%z
    3630                 particles(n)%speed_z = -particles(n)%speed_z
    3631                 IF ( use_sgs_for_particles  .AND. &
    3632                      particles(n)%rvar3 > 0.0 )  THEN
    3633                    particles(n)%rvar3 = -particles(n)%rvar3
    3634                 ENDIF
    3635                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3636                    particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
    3637                                                particle_tail_coordinates(1,3,nn)
    3638                 ENDIF
    3639              ENDIF
    3640           ENDIF
    3641           IF ( particles(n)%z < zw(0)  .AND.  particle_mask(n) )  THEN
    3642              IF ( ibc_par_b == 1 )  THEN
    3643 !
    3644 !--             Particle absorption
    3645                 particle_mask(n)  = .FALSE.
    3646                 deleted_particles = deleted_particles + 1
    3647                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3648                    tail_mask(nn) = .FALSE.
    3649                    deleted_tails = deleted_tails + 1
    3650                 ENDIF
    3651              ELSEIF ( ibc_par_b == 2 )  THEN
    3652 !
    3653 !--             Particle reflection
    3654                 particles(n)%z       = 2.0 * zw(0) - particles(n)%z
    3655                 particles(n)%speed_z = -particles(n)%speed_z
    3656                 IF ( use_sgs_for_particles  .AND. &
    3657                      particles(n)%rvar3 < 0.0 )  THEN
    3658                    particles(n)%rvar3 = -particles(n)%rvar3
    3659                 ENDIF
    3660                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3661                    particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
    3662                                                particle_tail_coordinates(1,3,nn)
    3663                 ENDIF
    3664                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    3665                    particle_tail_coordinates(1,3,nn) = 2.0 * zw(0) - &
    3666                                                particle_tail_coordinates(1,3,nn)
    3667                 ENDIF
    3668              ENDIF
    3669           ENDIF
    3670        ENDDO
    3671 
    3672 !    WRITE ( 9, * ) '*** advec_particles: ##7'
    3673 !    CALL local_flush( 9 )
    3674 !    nd = 0
    3675 !    DO  n = 1, number_of_particles
    3676 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    3677 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    3678 !       THEN
    3679 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3680 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    3681 !          CALL local_flush( 9 )
    3682 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3683 !       ENDIF
    3684 !    ENDDO
    3685 !    IF ( nd /= deleted_particles ) THEN
    3686 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    3687 !       CALL local_flush( 9 )
    3688 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    3689 !    ENDIF
     228       CALL lpm_boundary_conds( 'bottom/top' )
    3690229
    3691230!
     
    3693232!--    determine new number of particles
    3694233       IF ( number_of_particles > 0  .AND.  deleted_particles > 0 )  THEN
    3695           nn = 0
    3696           nd = 0
    3697           DO  n = 1, number_of_particles
    3698              IF ( particle_mask(n) )  THEN
    3699                 nn = nn + 1
    3700                 particles(nn) = particles(n)
    3701              ELSE
    3702                 nd = nd + 1
    3703              ENDIF
    3704           ENDDO
    3705 !       IF ( nd /= deleted_particles ) THEN
    3706 !          WRITE (9,*) '*** advec_part nd=',nd,' deleted_particles=',deleted_particles
    3707 !          CALL local_flush( 9 )
    3708 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3709 !       ENDIF
    3710 
    3711           number_of_particles = number_of_particles - deleted_particles
    3712 !
    3713 !--       Pack the tails, store the new tail ids and re-assign it to the
    3714 !--       respective
    3715 !--       particles
    3716           IF ( use_particle_tails )  THEN
    3717              nn = 0
    3718              nd = 0
    3719              DO  n = 1, number_of_tails
    3720                 IF ( tail_mask(n) )  THEN
    3721                    nn = nn + 1
    3722                    particle_tail_coordinates(:,:,nn) = &
    3723                                                 particle_tail_coordinates(:,:,n)
    3724                    new_tail_id(n) = nn
    3725                 ELSE
    3726                    nd = nd + 1
    3727 !                WRITE (9,*) '+++ n=',n,' (of ',number_of_tails,' #oftails)'
    3728 !                WRITE (9,*) '    id=',new_tail_id(n)
    3729 !                CALL local_flush( 9 )
    3730                 ENDIF
    3731              ENDDO
    3732           ENDIF
    3733 
    3734 !       IF ( nd /= deleted_tails  .AND.  use_particle_tails ) THEN
    3735 !          WRITE (9,*) '*** advec_part nd=',nd,' deleted_tails=',deleted_tails
    3736 !          CALL local_flush( 9 )
    3737 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3738 !       ENDIF
    3739 
    3740           number_of_tails = number_of_tails - deleted_tails
    3741 
    3742 !     nn = 0
    3743           DO  n = 1, number_of_particles
    3744              IF ( particles(n)%tail_id /= 0 )  THEN
    3745 !     nn = nn + 1
    3746 !     IF (  particles(n)%tail_id<0 .OR. particles(n)%tail_id > number_of_tails )  THEN
    3747 !        WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3748 !        WRITE (9,*) '    tail_id=',particles(n)%tail_id
    3749 !        WRITE (9,*) '    new_tail_id=', new_tail_id(particles(n)%tail_id), &
    3750 !                         ' of (',number_of_tails,')'
    3751 !        CALL local_flush( 9 )
    3752 !     ENDIF
    3753                 particles(n)%tail_id = new_tail_id(particles(n)%tail_id)
    3754              ENDIF
    3755           ENDDO
    3756 
    3757 !     IF ( nn /= number_of_tails  .AND.  use_particle_tails ) THEN
    3758 !        WRITE (9,*) '*** advec_part #of_tails=',number_of_tails,' nn=',nn
    3759 !        CALL local_flush( 9 )
    3760 !        DO  n = 1, number_of_particles
    3761 !           WRITE (9,*) 'prt# ',n,' tail_id=',particles(n)%tail_id, &
    3762 !                       ' x=',particles(n)%x, ' y=',particles(n)%y, &
    3763 !                       ' z=',particles(n)%z
    3764 !        ENDDO
    3765 !        CALL MPI_ABORT( comm2d, 9999, ierr )
    3766 !     ENDIF
    3767 
     234          CALL lpm_pack_arrays
    3768235       ENDIF
    3769236
    3770 !    IF ( number_of_particles /= number_of_tails )  THEN
    3771 !       WRITE (9,*) '--- advec_particles: #7'
    3772 !       WRITE (9,*) '    #of p=',number_of_particles,' #of t=',number_of_tails
    3773 !       CALL local_flush( 9 )
    3774 !    ENDIF
    3775 !    WRITE ( 9, * ) '*** advec_particles: ##8'
    3776 !    CALL local_flush( 9 )
    3777 !    DO  n = 1, number_of_particles
    3778 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    3779 !       THEN
    3780 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3781 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    3782 !          CALL local_flush( 9 )
    3783 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3784 !       ENDIF
    3785 !    ENDDO
    3786 
    3787 !    WRITE ( 9, * ) '*** advec_particles: ##9'
    3788 !    CALL local_flush( 9 )
    3789 !    DO  n = 1, number_of_particles
    3790 !       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
    3791 !       THEN
    3792 !          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
    3793 !          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
    3794 !          CALL local_flush( 9 )
    3795 !          CALL MPI_ABORT( comm2d, 9999, ierr )
    3796 !       ENDIF
    3797 !    ENDDO
    3798 
    3799 !
    3800 !--    Accumulate the number of particles transferred between the subdomains
    3801 #if defined( __parallel )
    3802        trlp_count_sum      = trlp_count_sum      + trlp_count
    3803        trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
    3804        trrp_count_sum      = trrp_count_sum      + trrp_count
    3805        trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
    3806        trsp_count_sum      = trsp_count_sum      + trsp_count
    3807        trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
    3808        trnp_count_sum      = trnp_count_sum      + trnp_count
    3809        trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
    3810 #endif
    3811 
    3812237       IF ( dt_3d_reached )  EXIT
    3813238
    3814 !
    3815 !--    Initialize variables for the next (sub-) timestep, i.e. for marking those
    3816 !--    particles to be deleted after the timestep
    3817        particle_mask     = .TRUE.
    3818        deleted_particles = 0
    3819        trlp_count_recv   = 0
    3820        trnp_count_recv   = 0
    3821        trrp_count_recv   = 0
    3822        trsp_count_recv   = 0
    3823        trlpt_count_recv  = 0
    3824        trnpt_count_recv  = 0
    3825        trrpt_count_recv  = 0
    3826        trspt_count_recv  = 0
    3827        IF ( use_particle_tails )  THEN
    3828           tail_mask     = .TRUE.
    3829        ENDIF
    3830        deleted_tails = 0
    3831 
    3832239    ENDDO   ! timestep loop
     240
    3833241
    3834242!
     
    3836244    time_sort_particles = time_sort_particles + dt_3d
    3837245    IF ( time_sort_particles >= dt_sort_particles )  THEN
    3838        CALL sort_particles
     246       CALL lpm_sort_arrays
    3839247       time_sort_particles = MOD( time_sort_particles, &
    3840248                                  MAX( dt_sort_particles, dt_3d ) )
    3841249    ENDIF
    3842250
    3843     IF ( cloud_droplets )  THEN
    3844 
    3845        CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'start' )
    3846 
    3847        ql = 0.0;  ql_v = 0.0;  ql_vp = 0.0
    3848 
    3849 !
    3850 !--    Calculate the liquid water content
    3851        DO  i = nxl, nxr
    3852           DO  j = nys, nyn
    3853              DO  k = nzb, nzt+1
    3854 
    3855 !
    3856 !--             Calculate the total volume in the boxes (ql_v, weighting factor
    3857 !--             has to beincluded)
    3858                 psi = prt_start_index(k,j,i)
    3859                 DO  n = psi, psi+prt_count(k,j,i)-1
    3860                    ql_v(k,j,i)  = ql_v(k,j,i)  + particles(n)%weight_factor *  &
    3861                                                  particles(n)%radius**3
    3862                 ENDDO
    3863 
    3864 !
    3865 !--             Calculate the liquid water content
    3866                 IF ( ql_v(k,j,i) /= 0.0 )  THEN
    3867                    ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi *           &
    3868                                            ql_v(k,j,i) /                       &
    3869                                            ( rho_surface * dx * dy * dz )
    3870 
    3871                    IF ( ql(k,j,i) < 0.0 ) THEN
    3872                       WRITE( message_string, * )  'LWC out of range: ' , &
    3873                                                   ql(k,j,i)
    3874                          CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 )
    3875                    ENDIF
    3876 
    3877                 ELSE
    3878                    ql(k,j,i) = 0.0
    3879                 ENDIF
    3880 
    3881              ENDDO
    3882           ENDDO
    3883        ENDDO
    3884 
    3885        CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'stop' )
    3886 
    3887     ENDIF
     251
     252!
     253!-- Calculate the new liquid water content for each grid box
     254    IF ( cloud_droplets )  CALL lpm_calc_liquid_water_content
     255
    3888256
    3889257!
     
    3892260!-- particle attribute (class) is then used for storing the particle radius
    3893261!-- class.
    3894     IF ( collision_kernel == 'none' )  CALL set_particle_attributes
     262    IF ( collision_kernel == 'none' )  CALL lpm_set_attributes
     263
    3895264
    3896265!
    3897266!-- Set particle attributes defined by the user
    3898     CALL user_particle_attributes
    3899 
    3900 !
    3901 !-- If necessary, add the actual particle positions to the particle tails
    3902     IF ( use_particle_tails )  THEN
    3903 
    3904        distance = 0.0
    3905        DO  n = 1, number_of_particles
    3906 
    3907           nn = particles(n)%tail_id
    3908 
    3909           IF ( nn /= 0 )  THEN
    3910 !
    3911 !--          Calculate the distance between the actual particle position and the
    3912 !--          next tailpoint
    3913 !             WRITE ( 9, * ) '*** advec_particles: ##10.1  nn=',nn
    3914 !             CALL local_flush( 9 )
    3915              IF ( minimum_tailpoint_distance /= 0.0 )  THEN
    3916                 distance = ( particle_tail_coordinates(1,1,nn) -      &
    3917                              particle_tail_coordinates(2,1,nn) )**2 + &
    3918                            ( particle_tail_coordinates(1,2,nn) -      &
    3919                              particle_tail_coordinates(2,2,nn) )**2 + &
    3920                            ( particle_tail_coordinates(1,3,nn) -      &
    3921                              particle_tail_coordinates(2,3,nn) )**2
    3922              ENDIF
    3923 !             WRITE ( 9, * ) '*** advec_particles: ##10.2'
    3924 !             CALL local_flush( 9 )
    3925 !
    3926 !--          First, increase the index of all existings tailpoints by one
    3927              IF ( distance >= minimum_tailpoint_distance )  THEN
    3928                 DO i = particles(n)%tailpoints, 1, -1
    3929                    particle_tail_coordinates(i+1,:,nn) = &
    3930                                                particle_tail_coordinates(i,:,nn)
    3931                 ENDDO
    3932 !
    3933 !--             Increase the counter which contains the number of tailpoints.
    3934 !--             This must always be smaller than the given maximum number of
    3935 !--             tailpoints because otherwise the index bounds of
    3936 !--             particle_tail_coordinates would be exceeded
    3937                 IF ( particles(n)%tailpoints < maximum_number_of_tailpoints-1 )&
    3938                 THEN
    3939                    particles(n)%tailpoints = particles(n)%tailpoints + 1
    3940                 ENDIF
    3941              ENDIF
    3942 !             WRITE ( 9, * ) '*** advec_particles: ##10.3'
    3943 !             CALL local_flush( 9 )
    3944 !
    3945 !--          In any case, store the new point at the beginning of the tail
    3946              particle_tail_coordinates(1,1,nn) = particles(n)%x
    3947              particle_tail_coordinates(1,2,nn) = particles(n)%y
    3948              particle_tail_coordinates(1,3,nn) = particles(n)%z
    3949              particle_tail_coordinates(1,4,nn) = particles(n)%class
    3950 !             WRITE ( 9, * ) '*** advec_particles: ##10.4'
    3951 !             CALL local_flush( 9 )
    3952 !
    3953 !--          Increase the age of the tailpoints
    3954              IF ( minimum_tailpoint_distance /= 0.0 )  THEN
    3955                 particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) =    &
    3956                    particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) + &
    3957                    dt_3d
    3958 !
    3959 !--             Delete the last tailpoint, if it has exceeded its maximum age
    3960                 IF ( particle_tail_coordinates(particles(n)%tailpoints,5,nn) > &
    3961                      maximum_tailpoint_age )  THEN
    3962                    particles(n)%tailpoints = particles(n)%tailpoints - 1
    3963                 ENDIF
    3964              ENDIF
    3965 !             WRITE ( 9, * ) '*** advec_particles: ##10.5'
    3966 !             CALL local_flush( 9 )
    3967 
    3968           ENDIF
    3969 
    3970        ENDDO
    3971 
    3972     ENDIF
    3973 !    WRITE ( 9, * ) '*** advec_particles: ##11'
    3974 !    CALL local_flush( 9 )
    3975 
    3976 !
    3977 !-- Write particle statistics on file
    3978     IF ( write_particle_statistics )  THEN
    3979        CALL check_open( 80 )
    3980 #if defined( __parallel )
    3981        WRITE ( 80, 8000 )  current_timestep_number+1, simulated_time+dt_3d, &
    3982                            number_of_particles, pleft, trlp_count_sum,      &
    3983                            trlp_count_recv_sum, pright, trrp_count_sum,     &
    3984                            trrp_count_recv_sum, psouth, trsp_count_sum,     &
    3985                            trsp_count_recv_sum, pnorth, trnp_count_sum,     &
    3986                            trnp_count_recv_sum, maximum_number_of_particles
    3987        CALL close_file( 80 )
    3988 #else
    3989        WRITE ( 80, 8000 )  current_timestep_number+1, simulated_time+dt_3d, &
    3990                            number_of_particles, maximum_number_of_particles
    3991 #endif
    3992     ENDIF
    3993 
    3994     CALL cpu_log( log_point(25), 'advec_particles', 'stop' )
    3995 
    3996 !
    3997 !-- Formats
    3998 8000 FORMAT (I6,1X,F7.2,4X,I6,5X,4(I3,1X,I4,'/',I4,2X),6X,I6)
    3999 
    4000  END SUBROUTINE advec_particles
    4001 
    4002 
    4003  SUBROUTINE allocate_prt_memory( number_of_new_particles )
    4004 
    4005 !------------------------------------------------------------------------------!
    4006 ! Description:
    4007 ! ------------
    4008 ! Extend particle memory
    4009 !------------------------------------------------------------------------------!
    4010 
    4011     USE particle_attributes
    4012 
    4013     IMPLICIT NONE
    4014 
    4015     INTEGER ::  new_maximum_number, number_of_new_particles
    4016 
    4017     LOGICAL, DIMENSION(:), ALLOCATABLE ::  tmp_particle_mask
    4018 
    4019     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  tmp_particles
    4020 
    4021 
    4022     new_maximum_number = maximum_number_of_particles + &
    4023                    MAX( 5*number_of_new_particles, number_of_initial_particles )
    4024 
    4025     IF ( write_particle_statistics )  THEN
    4026        CALL check_open( 80 )
    4027        WRITE ( 80, '(''*** Request: '', I7, '' new_maximum_number(prt)'')' ) &
    4028                             new_maximum_number
    4029        CALL close_file( 80 )
    4030     ENDIF
    4031 
    4032     ALLOCATE( tmp_particles(maximum_number_of_particles), &
    4033               tmp_particle_mask(maximum_number_of_particles) )
    4034     tmp_particles     = particles
    4035     tmp_particle_mask = particle_mask
    4036 
    4037     DEALLOCATE( particles, particle_mask )
    4038     ALLOCATE( particles(new_maximum_number), &
    4039               particle_mask(new_maximum_number) )
    4040     maximum_number_of_particles = new_maximum_number
    4041 
    4042     particles(1:number_of_particles) = tmp_particles(1:number_of_particles)
    4043     particle_mask(1:number_of_particles) = &
    4044                                   tmp_particle_mask(1:number_of_particles)
    4045     particle_mask(number_of_particles+1:maximum_number_of_particles) = .TRUE.
    4046     DEALLOCATE( tmp_particles, tmp_particle_mask )
    4047 
    4048  END SUBROUTINE allocate_prt_memory
    4049 
    4050 
    4051  SUBROUTINE allocate_tail_memory( number_of_new_tails )
    4052 
    4053 !------------------------------------------------------------------------------!
    4054 ! Description:
    4055 ! ------------
    4056 ! Extend tail memory
    4057 !------------------------------------------------------------------------------!
    4058 
    4059     USE particle_attributes
    4060 
    4061     IMPLICIT NONE
    4062 
    4063     INTEGER ::  new_maximum_number, number_of_new_tails
    4064 
    4065     LOGICAL, DIMENSION(maximum_number_of_tails) ::  tmp_tail_mask
    4066 
    4067     REAL, DIMENSION(maximum_number_of_tailpoints,5,maximum_number_of_tails) :: &
    4068                                                     tmp_tail
    4069 
    4070 
    4071     new_maximum_number = maximum_number_of_tails + &
    4072                          MAX( 5*number_of_new_tails, number_of_initial_tails )
    4073 
    4074     IF ( write_particle_statistics )  THEN
    4075        CALL check_open( 80 )
    4076        WRITE ( 80, '(''*** Request: '', I5, '' new_maximum_number(tails)'')' ) &
    4077                             new_maximum_number
    4078        CALL close_file( 80 )
    4079     ENDIF
    4080     WRITE (9,*) '*** Request: ',new_maximum_number,' new_maximum_number(tails)'
    4081 !    CALL local_flush( 9 )
    4082 
    4083     tmp_tail(:,:,1:number_of_tails)  = &
    4084                                 particle_tail_coordinates(:,:,1:number_of_tails)
    4085     tmp_tail_mask(1:number_of_tails) = tail_mask(1:number_of_tails)
    4086 
    4087     DEALLOCATE( new_tail_id, particle_tail_coordinates, tail_mask )
    4088     ALLOCATE( new_tail_id(new_maximum_number),                          &
    4089               particle_tail_coordinates(maximum_number_of_tailpoints,5, &
    4090               new_maximum_number),                                      &
    4091               tail_mask(new_maximum_number) )
    4092     maximum_number_of_tails = new_maximum_number
    4093 
    4094     particle_tail_coordinates = 0.0
    4095     particle_tail_coordinates(:,:,1:number_of_tails) = &
    4096                                                  tmp_tail(:,:,1:number_of_tails)
    4097     tail_mask(1:number_of_tails) = tmp_tail_mask(1:number_of_tails)
    4098     tail_mask(number_of_tails+1:maximum_number_of_tails) = .TRUE.
    4099 
    4100  END SUBROUTINE allocate_tail_memory
    4101 
    4102 
    4103  SUBROUTINE output_particles_netcdf
    4104 #if defined( __netcdf )
    4105 
    4106     USE control_parameters
    4107     USE netcdf_control
    4108     USE particle_attributes
    4109 
    4110     IMPLICIT NONE
    4111 
    4112 
    4113     CALL check_open( 108 )
    4114 
    4115 !
    4116 !-- Update the NetCDF time axis
    4117     prt_time_count = prt_time_count + 1
    4118 
    4119     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_time_prt, (/ simulated_time /), &
    4120                             start = (/ prt_time_count /), count = (/ 1 /) )
    4121     CALL handle_netcdf_error( 'output_particles_netcdf', 1 )
    4122 
    4123 !
    4124 !-- Output the real number of particles used
    4125     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_rnop_prt, &
    4126                             (/ number_of_particles /),   &
    4127                             start = (/ prt_time_count /), count = (/ 1 /) )
    4128     CALL handle_netcdf_error( 'output_particles_netcdf', 2 )
    4129 
    4130 !
    4131 !-- Output all particle attributes
    4132     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(1), particles%age,         &
    4133                             start = (/ 1, prt_time_count /),                  &
    4134                             count = (/ maximum_number_of_particles /) )
    4135     CALL handle_netcdf_error( 'output_particles_netcdf', 3 )
    4136 
    4137     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%dvrp_psize,  &
    4138                             start = (/ 1, prt_time_count /),                  &
    4139                             count = (/ maximum_number_of_particles /) )
    4140     CALL handle_netcdf_error( 'output_particles_netcdf', 4 )
    4141 
    4142     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(3), particles%origin_x,    &
    4143                             start = (/ 1, prt_time_count /),                  &
    4144                             count = (/ maximum_number_of_particles /) )
    4145     CALL handle_netcdf_error( 'output_particles_netcdf', 5 )
    4146 
    4147     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(4), particles%origin_y,    &
    4148                             start = (/ 1, prt_time_count /),                  &
    4149                             count = (/ maximum_number_of_particles /) )
    4150     CALL handle_netcdf_error( 'output_particles_netcdf', 6 )
    4151 
    4152     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(5), particles%origin_z,    &
    4153                             start = (/ 1, prt_time_count /),                  &
    4154                             count = (/ maximum_number_of_particles /) )
    4155     CALL handle_netcdf_error( 'output_particles_netcdf', 7 )
    4156 
    4157     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(6), particles%radius,      &
    4158                             start = (/ 1, prt_time_count /),                  &
    4159                             count = (/ maximum_number_of_particles /) )
    4160     CALL handle_netcdf_error( 'output_particles_netcdf', 8 )
    4161 
    4162     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(7), particles%speed_x,     &
    4163                             start = (/ 1, prt_time_count /),                  &
    4164                             count = (/ maximum_number_of_particles /) )
    4165     CALL handle_netcdf_error( 'output_particles_netcdf', 9 )
    4166 
    4167     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(8), particles%speed_y,     &
    4168                             start = (/ 1, prt_time_count /),                  &
    4169                             count = (/ maximum_number_of_particles /) )
    4170     CALL handle_netcdf_error( 'output_particles_netcdf', 10 )
    4171 
    4172     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(9), particles%speed_z,     &
    4173                             start = (/ 1, prt_time_count /),                  &
    4174                             count = (/ maximum_number_of_particles /) )
    4175     CALL handle_netcdf_error( 'output_particles_netcdf', 11 )
    4176 
    4177     nc_stat = NF90_PUT_VAR( id_set_prt,id_var_prt(10),particles%weight_factor,&
    4178                             start = (/ 1, prt_time_count /),                  &
    4179                             count = (/ maximum_number_of_particles /) )
    4180     CALL handle_netcdf_error( 'output_particles_netcdf', 12 )
    4181 
    4182     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(11), particles%x,          &
    4183                             start = (/ 1, prt_time_count /),                  &
    4184                             count = (/ maximum_number_of_particles /) )
    4185     CALL handle_netcdf_error( 'output_particles_netcdf', 13 )
    4186 
    4187     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(12), particles%y,          &
    4188                             start = (/ 1, prt_time_count /),                  &
    4189                             count = (/ maximum_number_of_particles /) )
    4190     CALL handle_netcdf_error( 'output_particles_netcdf', 14 )
    4191 
    4192     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(13), particles%z,          &
    4193                             start = (/ 1, prt_time_count /),                  &
    4194                             count = (/ maximum_number_of_particles /) )
    4195     CALL handle_netcdf_error( 'output_particles_netcdf', 15 )
    4196 
    4197     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%class,      &
    4198                             start = (/ 1, prt_time_count /),                  &
    4199                             count = (/ maximum_number_of_particles /) )
    4200     CALL handle_netcdf_error( 'output_particles_netcdf', 16 )
    4201 
    4202     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(15), particles%group,      &
    4203                             start = (/ 1, prt_time_count /),                  &
    4204                             count = (/ maximum_number_of_particles /) )
    4205     CALL handle_netcdf_error( 'output_particles_netcdf', 17 )
    4206 
    4207     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16), particles%tailpoints, &
    4208                             start = (/ 1, prt_time_count /),                  &
    4209                             count = (/ maximum_number_of_particles /) )
    4210     CALL handle_netcdf_error( 'output_particles_netcdf', 18 )
    4211 
    4212     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%tail_id,    &
    4213                             start = (/ 1, prt_time_count /),                  &
    4214                             count = (/ maximum_number_of_particles /) )
    4215     CALL handle_netcdf_error( 'output_particles_netcdf', 19 )
    4216 
    4217 #endif
    4218  END SUBROUTINE output_particles_netcdf
    4219 
    4220 
    4221  SUBROUTINE write_particles
    4222 
    4223 !------------------------------------------------------------------------------!
    4224 ! Description:
    4225 ! ------------
    4226 ! Write particle data on restart file
    4227 !------------------------------------------------------------------------------!
    4228 
    4229     USE control_parameters
    4230     USE particle_attributes
    4231     USE pegrid
    4232 
    4233     IMPLICIT NONE
    4234 
    4235     CHARACTER (LEN=10) ::  particle_binary_version
    4236     INTEGER ::  i
    4237 
    4238 !
    4239 !-- First open the output unit.
    4240     IF ( myid_char == '' )  THEN
    4241        OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT'//myid_char, &
    4242                   FORM='UNFORMATTED')
    4243     ELSE
    4244        IF ( myid == 0 )  CALL local_system( 'mkdir PARTICLE_RESTART_DATA_OUT' )
    4245 #if defined( __parallel )
    4246 !
    4247 !--    Set a barrier in order to allow that thereafter all other processors
    4248 !--    in the directory created by PE0 can open their file
    4249        CALL MPI_BARRIER( comm2d, ierr )
    4250 #endif
    4251        OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT/'//myid_char, &
    4252                   FORM='UNFORMATTED' )
    4253     ENDIF
    4254 
    4255     DO  i = 0, io_blocks-1
    4256 
    4257        IF ( i == io_group )  THEN
    4258 
    4259 !
    4260 !--       Write the version number of the binary format.
    4261 !--       Attention: After changes to the following output commands the version
    4262 !--       ---------  number of the variable particle_binary_version must be
    4263 !--                  changed! Also, the version number and the list of arrays
    4264 !--                  to be read in init_particles must be adjusted accordingly.
    4265           particle_binary_version = '3.0'
    4266           WRITE ( 90 )  particle_binary_version
    4267 
    4268 !
    4269 !--       Write some particle parameters, the size of the particle arrays as
    4270 !--       well as other dvrp-plot variables.
    4271           WRITE ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,              &
    4272                         maximum_number_of_particles,                           &
    4273                         maximum_number_of_tailpoints, maximum_number_of_tails, &
    4274                         number_of_initial_particles, number_of_particles,      &
    4275                         number_of_particle_groups, number_of_tails,            &
    4276                         particle_groups, time_prel, time_write_particle_data,  &
    4277                         uniform_particles
    4278 
    4279           IF ( number_of_initial_particles /= 0 ) WRITE ( 90 ) initial_particles
    4280 
    4281           WRITE ( 90 )  prt_count, prt_start_index
    4282           WRITE ( 90 )  particles
    4283 
    4284           IF ( use_particle_tails )  THEN
    4285              WRITE ( 90 )  particle_tail_coordinates
    4286           ENDIF
    4287 
    4288           CLOSE ( 90 )
    4289 
    4290        ENDIF
    4291 
    4292 #if defined( __parallel )
    4293        CALL MPI_BARRIER( comm2d, ierr )
    4294 #endif
    4295 
    4296     ENDDO
    4297 
    4298  END SUBROUTINE write_particles
    4299 
    4300 
    4301 
    4302  SUBROUTINE collision_efficiency( mean_r, r, e)
    4303 !------------------------------------------------------------------------------!
    4304 ! Description:
    4305 ! ------------
    4306 ! Interpolate collision efficiency from table
    4307 !------------------------------------------------------------------------------!
    4308 
    4309     IMPLICIT NONE
    4310 
    4311     INTEGER       ::  i, j, k
    4312 
    4313     LOGICAL, SAVE ::  first = .TRUE.
    4314 
    4315     REAL          ::  aa, bb, cc, dd, dx, dy, e, gg, mean_r, mean_rm, r, rm, &
    4316                       x, y
    4317 
    4318     REAL, DIMENSION(1:9), SAVE      ::  collected_r = 0.0
    4319     REAL, DIMENSION(1:19), SAVE     ::  collector_r = 0.0
    4320     REAL, DIMENSION(1:9,1:19), SAVE ::  ef = 0.0
    4321 
    4322     mean_rm = mean_r * 1.0E06
    4323     rm      = r      * 1.0E06
    4324 
    4325     IF ( first )  THEN
    4326        collected_r = (/ 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25.0 /)
    4327        collector_r = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0, 150.0,&
    4328                         200.0, 300.0, 400.0, 500.0, 600.0, 1000.0, 1400.0,     &
    4329                         1800.0, 2400.0, 3000.0 /)
    4330        ef(:,1) = (/0.017, 0.027, 0.037, 0.052, 0.052, 0.052, 0.052, 0.0,  0.0 /)
    4331        ef(:,2) = (/0.001, 0.016, 0.027, 0.060, 0.12,  0.17,  0.17,  0.17, 0.0 /)
    4332        ef(:,3) = (/0.001, 0.001, 0.02,  0.13,  0.28,  0.37,  0.54,  0.55, 0.47/)
    4333        ef(:,4) = (/0.001, 0.001, 0.02,  0.23,  0.4,   0.55,  0.7,   0.75, 0.75/)
    4334        ef(:,5) = (/0.01,  0.01,  0.03,  0.3,   0.4,   0.58,  0.73,  0.75, 0.79/)
    4335        ef(:,6) = (/0.01,  0.01,  0.13,  0.38,  0.57,  0.68,  0.80,  0.86, 0.91/)
    4336        ef(:,7) = (/0.01,  0.085, 0.23,  0.52,  0.68,  0.76,  0.86,  0.92, 0.95/)
    4337        ef(:,8) = (/0.01,  0.14,  0.32,  0.60,  0.73,  0.81,  0.90,  0.94, 0.96/)
    4338        ef(:,9) = (/0.025, 0.25,  0.43,  0.66,  0.78,  0.83,  0.92,  0.95, 0.96/)
    4339        ef(:,10)= (/0.039, 0.3,   0.46,  0.69,  0.81,  0.87,  0.93,  0.95, 0.96/)
    4340        ef(:,11)= (/0.095, 0.33,  0.51,  0.72,  0.82,  0.87,  0.93,  0.96, 0.97/)
    4341        ef(:,12)= (/0.098, 0.36,  0.51,  0.73,  0.83,  0.88,  0.93,  0.96, 0.97/)
    4342        ef(:,13)= (/0.1,   0.36,  0.52,  0.74,  0.83,  0.88,  0.93,  0.96, 0.97/)
    4343        ef(:,14)= (/0.17,  0.4,   0.54,  0.72,  0.83,  0.88,  0.94,  0.98, 1.0 /)
    4344        ef(:,15)= (/0.15,  0.37,  0.52,  0.74,  0.82,  0.88,  0.94,  0.98, 1.0 /)
    4345        ef(:,16)= (/0.11,  0.34,  0.49,  0.71,  0.83,  0.88,  0.94,  0.95, 1.0 /)
    4346        ef(:,17)= (/0.08,  0.29,  0.45,  0.68,  0.8,   0.86,  0.96,  0.94, 1.0 /)
    4347        ef(:,18)= (/0.04,  0.22,  0.39,  0.62,  0.75,  0.83,  0.92,  0.96, 1.0 /)
    4348        ef(:,19)= (/0.02,  0.16,  0.33,  0.55,  0.71,  0.81,  0.90,  0.94, 1.0 /)
    4349     ENDIF
    4350 
    4351     DO  k = 1, 8
    4352        IF ( collected_r(k) <= mean_rm )  i = k
    4353     ENDDO
    4354 
    4355     DO  k = 1, 18
    4356        IF ( collector_r(k) <= rm )  j = k
    4357     ENDDO
    4358 
    4359     IF ( rm < 10.0 )  THEN
    4360        e = 0.0
    4361     ELSEIF ( mean_rm < 2.0 )  THEN
    4362        e = 0.001
    4363     ELSEIF ( mean_rm >= 25.0 )  THEN
    4364        IF( j <= 2 )  e = 0.0
    4365        IF( j == 3 )  e = 0.47
    4366        IF( j == 4 )  e = 0.8
    4367        IF( j == 5 )  e = 0.9
    4368        IF( j >=6  )  e = 1.0
    4369     ELSEIF ( rm >= 3000.0 )  THEN
    4370        IF( i == 1 )  e = 0.02
    4371        IF( i == 2 )  e = 0.16
    4372        IF( i == 3 )  e = 0.33
    4373        IF( i == 4 )  e = 0.55
    4374        IF( i == 5 )  e = 0.71
    4375        IF( i == 6 )  e = 0.81
    4376        IF( i == 7 )  e = 0.90
    4377        IF( i >= 8 )  e = 0.94
    4378     ELSE
    4379        x  = mean_rm - collected_r(i)
    4380        y  = rm - collector_r(j)
    4381        dx = collected_r(i+1) - collected_r(i)
    4382        dy = collector_r(j+1) - collector_r(j)
    4383        aa = x**2 + y**2
    4384        bb = ( dx - x )**2 + y**2
    4385        cc = x**2 + ( dy - y )**2
    4386        dd = ( dx - x )**2 + ( dy - y )**2
    4387        gg = aa + bb + cc + dd
    4388 
    4389        e = ( (gg-aa)*ef(i,j) + (gg-bb)*ef(i+1,j) + (gg-cc)*ef(i,j+1) + &
    4390              (gg-dd)*ef(i+1,j+1) ) / (3.0*gg)
    4391     ENDIF
    4392 
    4393  END SUBROUTINE collision_efficiency
    4394 
    4395 
    4396 
    4397  SUBROUTINE sort_particles
    4398 
    4399 !------------------------------------------------------------------------------!
    4400 ! Description:
    4401 ! ------------
    4402 ! Sort particles in the sequence the grid boxes are stored in memory
    4403 !------------------------------------------------------------------------------!
    4404 
    4405     USE arrays_3d
    4406     USE control_parameters
    4407     USE cpulog
    4408     USE grid_variables
    4409     USE indices
    4410     USE interfaces
    4411     USE particle_attributes
    4412 
    4413     IMPLICIT NONE
    4414 
    4415     INTEGER ::  i, ilow, j, k, n
    4416 
    4417     TYPE(particle_type), DIMENSION(:), POINTER ::  particles_temp
    4418 
    4419 
    4420     CALL cpu_log( log_point_s(47), 'sort_particles', 'start' )
    4421 
    4422 !
    4423 !-- Initialize counters and set pointer of the temporary array into which
    4424 !-- particles are sorted to free memory
    4425     prt_count  = 0
    4426     sort_count = sort_count +1
    4427 
    4428     SELECT CASE ( MOD( sort_count, 2 ) )
    4429 
    4430        CASE ( 0 )
    4431 
    4432           particles_temp => part_1
    4433 !          part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    4434 !                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    4435 !                                  0.0, 0, 0, 0, 0 )
    4436 
    4437        CASE ( 1 )
    4438 
    4439           particles_temp => part_2
    4440 !          part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    4441 !                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    4442 !                                  0.0, 0, 0, 0, 0 )
    4443 
    4444     END SELECT
    4445 
    4446 !
    4447 !-- Count the particles per gridbox
    4448     DO  n = 1, number_of_particles
    4449 
    4450        i = ( particles(n)%x + 0.5 * dx ) * ddx
    4451        j = ( particles(n)%y + 0.5 * dy ) * ddy
    4452        k = particles(n)%z / dz + 1 + offset_ocean_nzt
    4453            ! only exact if equidistant
    4454 
    4455        prt_count(k,j,i) = prt_count(k,j,i) + 1
    4456 
    4457        IF ( i < nxl .OR. i > nxr .OR. j < nys .OR. j > nyn .OR. k < nzb+1 .OR. &
    4458             k > nzt )  THEN
    4459           WRITE( message_string, * ) ' particle out of range: i=', i, ' j=', &
    4460                           j, ' k=', k,                                       &
    4461                           ' nxl=', nxl, ' nxr=', nxr,                        &
    4462                           ' nys=', nys, ' nyn=', nyn,                        &
    4463                           ' nzb=', nzb, ' nzt=', nzt
    4464          CALL message( 'sort_particles', 'PA0149', 1, 2, 0, 6, 0 )
    4465        ENDIF
    4466 
    4467     ENDDO
    4468 
    4469 !
    4470 !-- Calculate the lower indices of those ranges of the particles-array
    4471 !-- containing particles which belong to the same gridpox i,j,k
    4472     ilow = 1
    4473     DO  i = nxl, nxr
    4474        DO  j = nys, nyn
    4475           DO  k = nzb+1, nzt
    4476              prt_start_index(k,j,i) = ilow
    4477              ilow = ilow + prt_count(k,j,i)
    4478           ENDDO
    4479        ENDDO
    4480     ENDDO
    4481 
    4482 !
    4483 !-- Sorting the particles
    4484     DO  n = 1, number_of_particles
    4485 
    4486        i = ( particles(n)%x + 0.5 * dx ) * ddx
    4487        j = ( particles(n)%y + 0.5 * dy ) * ddy
    4488        k = particles(n)%z / dz + 1 + offset_ocean_nzt
    4489            ! only exact if equidistant
    4490 
    4491        particles_temp(prt_start_index(k,j,i)) = particles(n)
    4492 
    4493        prt_start_index(k,j,i) = prt_start_index(k,j,i) + 1
    4494 
    4495     ENDDO
    4496 
    4497 !
    4498 !-- Redirect the particles pointer to the sorted array
    4499     SELECT CASE ( MOD( sort_count, 2 ) )
    4500 
    4501        CASE ( 0 )
    4502 
    4503           particles => part_1
    4504 
    4505        CASE ( 1 )
    4506 
    4507           particles => part_2
    4508 
    4509     END SELECT
    4510 
    4511 !
    4512 !-- Reset the index array to the actual start position
    4513     DO  i = nxl, nxr
    4514        DO  j = nys, nyn
    4515           DO  k = nzb+1, nzt
    4516              prt_start_index(k,j,i) = prt_start_index(k,j,i) - prt_count(k,j,i)
    4517           ENDDO
    4518        ENDDO
    4519     ENDDO
    4520 
    4521     CALL cpu_log( log_point_s(47), 'sort_particles', 'stop' )
    4522 
    4523  END SUBROUTINE sort_particles
     267    CALL user_lpm_set_attributes
     268
     269
     270!
     271!-- If required, add the current particle positions to the particle tails
     272    IF ( use_particle_tails )  CALL lpm_extend_tails
     273
     274
     275!
     276!-- Write particle statistics (in particular the number of particles
     277!-- exchanged between the subdomains) on file
     278    IF ( write_particle_statistics )  CALL lpm_write_exchange_statistics
     279
     280    CALL cpu_log( log_point(25), 'lpm', 'stop' )
     281
     282
     283 END SUBROUTINE lpm
  • palm/trunk/SOURCE/lpm_boundary_conds.f90

    r848 r849  
    1  SUBROUTINE particle_boundary_conds
     1 SUBROUTINE lpm_boundary_conds( range )
    22
    33!------------------------------------------------------------------------------!
    44! Current revisions:
    55! -----------------
    6 !
     6! routine renamed lpm_boundary_conds, bottom and top boundary conditions
     7! included
    78!
    89! Former revisions:
     
    2021! Description:
    2122! ------------
    22 ! Calculates the reflection of particles from vertical walls.
    23 ! Routine developed by Jin Zhang (2006-2007)
     23! Boundary conditions for the Lagrangian particles.
     24! The routine consists of two different parts. One handles the bottom (flat)
     25! and top boundary. In this part, also particles which exceeded their lifetime
     26! are deleted.
     27! The other part handles the reflection of particles from vertical walls.
     28! This part was developed by Jin Zhang during 2006-2007.
     29!
    2430! To do: Code structure for finding the t_index values and for checking the
    25 !        reflection conditions is basically the same for all four cases, so it
     31! -----  reflection conditions is basically the same for all four cases, so it
    2632!        should be possible to further simplify/shorten it.
    27 ! THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR! (see offset_ocean_*)
     33!
     34! THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
     35! (see offset_ocean_*)
    2836!------------------------------------------------------------------------------!
    2937
     38    USE arrays_3d
    3039    USE control_parameters
    3140    USE cpulog
     
    3847    IMPLICIT NONE
    3948
     49    CHARACTER (LEN=*) ::  range
     50
    4051    INTEGER ::  i, inc, ir, i1, i2, i3, i5, j, jr, j1, j2, j3, j5, k, k1, k2, &
    41                 k3, k5, n, t_index, t_index_number
     52                k3, k5, n, nn, t_index, t_index_number
    4253
    4354    LOGICAL ::  reflect_x, reflect_y, reflect_z
     
    4859    REAL ::  t(1:200)
    4960
    50     CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' )
    51 
    52 
    53 
    54     reflect_x = .FALSE.
    55     reflect_y = .FALSE.
    56     reflect_z = .FALSE.
    57 
    58     DO  n = 1, number_of_particles
    59 
    60        dt_particle = particles(n)%age - particles(n)%age_m
    61 
    62        i2 = ( particles(n)%x + 0.5 * dx ) * ddx
    63        j2 = ( particles(n)%y + 0.5 * dy ) * ddy
    64        k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
    65 
    66        prt_x = particles(n)%x
    67        prt_y = particles(n)%y
    68        prt_z = particles(n)%z
    69 
    70 !
    71 !--    If the particle position is below the surface, it has to be reflected
    72        IF ( k2 <= nzb_s_inner(j2,i2)  .AND.  nzb_s_inner(j2,i2) /=0 )  THEN
    73 
    74           pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
    75           pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
    76           pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
    77           i1 = ( pos_x_old + 0.5 * dx ) * ddx
    78           j1 = ( pos_y_old + 0.5 * dy ) * ddy
    79           k1 = pos_z_old / dz + offset_ocean_nzt_m1
    80 
    81 !
    82 !--       Case 1
    83           IF ( particles(n)%x > pos_x_old  .AND.  particles(n)%y > pos_y_old ) &
     61
     62
     63    IF ( range == 'bottom/top' )  THEN
     64
     65!
     66!--    Apply boundary conditions to those particles that have crossed the top or
     67!--    bottom boundary and delete those particles, which are older than allowed
     68       DO  n = 1, number_of_particles
     69
     70          nn = particles(n)%tail_id
     71
     72!
     73!--       Stop if particles have moved further than the length of one
     74!--       PE subdomain (newly released particles have age = age_m!)
     75          IF ( particles(n)%age /= particles(n)%age_m )  THEN
     76             IF ( ABS(particles(n)%speed_x) >                                  &
     77                  ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m)  .OR. &
     78                  ABS(particles(n)%speed_y) >                                  &
     79                  ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) )  THEN
     80
     81                  WRITE( message_string, * )  'particle too fast.  n = ',  n
     82                  CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
     83             ENDIF
     84          ENDIF
     85
     86          IF ( particles(n)%age > particle_maximum_age  .AND.  &
     87               particle_mask(n) )                              &
    8488          THEN
    85              t_index = 1
    86 
    87              DO  i = i1, i2
    88                 xline      = i * dx + 0.5 * dx
    89                 t(t_index) = ( xline - pos_x_old ) / &
    90                              ( particles(n)%x - pos_x_old )
    91                 t_index    = t_index + 1
    92              ENDDO
    93 
    94              DO  j = j1, j2
    95                 yline      = j * dy + 0.5 * dy
    96                 t(t_index) = ( yline - pos_y_old ) / &
    97                              ( particles(n)%y - pos_y_old )
    98                 t_index    = t_index + 1
    99              ENDDO
    100 
    101              IF ( particles(n)%z < pos_z_old )  THEN
    102                 DO  k = k1, k2 , -1
    103                    zline      = k * dz
    104                    t(t_index) = ( pos_z_old - zline ) / &
    105                                 ( pos_z_old - particles(n)%z )
    106                    t_index    = t_index + 1
    107                 ENDDO
    108              ENDIF
    109 
    110              t_index_number = t_index - 1
    111 
    112 !
    113 !--          Sorting t
    114              inc = 1
    115              jr  = 1
    116              DO WHILE ( inc <= t_index_number )
    117                 inc = 3 * inc + 1
    118              ENDDO
    119 
    120              DO WHILE ( inc > 1 )
    121                 inc = inc / 3
    122                 DO  ir = inc+1, t_index_number
    123                    tmp_t = t(ir)
    124                    jr    = ir
    125                    DO WHILE ( t(jr-inc) > tmp_t )
    126                       t(jr) = t(jr-inc)
    127                       jr    = jr - inc
    128                       IF ( jr <= inc )  EXIT
    129                    ENDDO
    130                    t(jr) = tmp_t
    131                 ENDDO
    132              ENDDO
    133 
    134       case1: DO  t_index = 1, t_index_number
    135 
    136                 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    137                 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    138                 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    139 
    140                 i3 = ( pos_x + 0.5 * dx ) * ddx   
    141                 j3 = ( pos_y + 0.5 * dy ) * ddy
    142                 k3 = pos_z / dz + offset_ocean_nzt_m1
    143 
    144                 i5 = pos_x * ddx
    145                 j5 = pos_y * ddy
    146                 k5 = pos_z / dz + offset_ocean_nzt_m1
    147 
    148                 IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    149                      nzb_s_inner(j5,i3) /= 0 )  THEN
    150 
    151                    IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    152                         k3 == nzb_s_inner(j5,i3) )  THEN
    153                       reflect_z = .TRUE.
    154                       EXIT case1
    155                    ENDIF
    156 
    157                 ENDIF
    158 
    159                 IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    160                      nzb_s_inner(j3,i5) /= 0 )  THEN
    161 
    162                    IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    163                         k3 == nzb_s_inner(j3,i5) )  THEN
    164                       reflect_z = .TRUE.
    165                       EXIT case1
    166                    ENDIF
    167 
    168                 ENDIF
    169 
    170                 IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    171                      nzb_s_inner(j3,i3) /= 0 )  THEN
    172 
    173                    IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    174                         k3 == nzb_s_inner(j3,i3) )  THEN
    175                       reflect_z = .TRUE.
    176                       EXIT case1
    177                    ENDIF
    178 
    179                    IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
    180                         pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    181                       reflect_y = .TRUE.
    182                       EXIT case1
    183                    ENDIF
    184 
    185                    IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
    186                         pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    187                       reflect_x = .TRUE.
    188                       EXIT case1
    189                    ENDIF
    190 
    191                 ENDIF
    192 
    193              ENDDO case1
    194 
    195 !
    196 !--       Case 2
    197           ELSEIF ( particles(n)%x > pos_x_old  .AND. &
    198                    particles(n)%y < pos_y_old )  THEN
    199 
    200              t_index = 1
    201 
    202              DO  i = i1, i2
    203                 xline      = i * dx + 0.5 * dx
    204                 t(t_index) = ( xline - pos_x_old ) / &
    205                              ( particles(n)%x - pos_x_old )
    206                 t_index    = t_index + 1
    207              ENDDO
    208 
    209              DO  j = j1, j2, -1
    210                 yline      = j * dy - 0.5 * dy
    211                 t(t_index) = ( pos_y_old - yline ) / &
    212                              ( pos_y_old - particles(n)%y )
    213                 t_index    = t_index + 1
    214              ENDDO
    215 
    216              IF ( particles(n)%z < pos_z_old )  THEN
    217                 DO  k = k1, k2 , -1
    218                    zline      = k * dz
    219                    t(t_index) = ( pos_z_old - zline ) / &
    220                                 ( pos_z_old - particles(n)%z )
    221                    t_index    = t_index + 1
    222                 ENDDO
    223              ENDIF
    224              t_index_number = t_index-1
    225 
    226 !
    227 !--          Sorting t
    228              inc = 1
    229              jr  = 1
    230              DO WHILE ( inc <= t_index_number )
    231                 inc = 3 * inc + 1
    232              ENDDO
    233 
    234              DO WHILE ( inc > 1 )
    235                 inc = inc / 3
    236                 DO  ir = inc+1, t_index_number
    237                    tmp_t = t(ir)
    238                    jr    = ir
    239                    DO WHILE ( t(jr-inc) > tmp_t )
    240                       t(jr) = t(jr-inc)
    241                       jr    = jr - inc
    242                       IF ( jr <= inc )  EXIT
    243                    ENDDO
    244                    t(jr) = tmp_t
    245                 ENDDO
    246              ENDDO
    247 
    248       case2: DO  t_index = 1, t_index_number
    249 
    250                 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    251                 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    252                 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    253 
    254                 i3 = ( pos_x + 0.5 * dx ) * ddx
    255                 j3 = ( pos_y + 0.5 * dy ) * ddy
    256                 k3 = pos_z / dz + offset_ocean_nzt_m1
    257 
    258                 i5 = pos_x * ddx
    259                 j5 = pos_y * ddy
    260                 k5 = pos_z / dz + offset_ocean_nzt_m1
    261 
    262                 IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    263                      nzb_s_inner(j3,i5) /= 0 )  THEN
    264 
    265                    IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    266                         k3 == nzb_s_inner(j3,i5) )  THEN
    267                       reflect_z = .TRUE.
    268                       EXIT case2
    269                    ENDIF
    270 
    271                 ENDIF
    272 
    273                 IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    274                      nzb_s_inner(j3,i3) /= 0 )  THEN
    275 
    276                    IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    277                         k3 == nzb_s_inner(j3,i3) )  THEN
    278                       reflect_z = .TRUE.
    279                       EXIT case2
    280                    ENDIF
    281 
    282                    IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
    283                         pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    284                       reflect_x = .TRUE.
    285                       EXIT case2
    286                    ENDIF
    287 
    288                 ENDIF
    289 
    290                 IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    291                      nzb_s_inner(j5,i3) /= 0 )  THEN
    292 
    293                    IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    294                         k3 == nzb_s_inner(j5,i3) )  THEN
    295                       reflect_z = .TRUE.
    296                       EXIT case2
    297                    ENDIF
    298 
    299                    IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
    300                         pos_z < nzb_s_inner(j5,i3) * dz )  THEN
    301                       reflect_y = .TRUE.
    302                       EXIT case2
    303                    ENDIF
    304 
    305                 ENDIF
    306 
    307              ENDDO case2
    308 
    309 !
    310 !--       Case 3
    311           ELSEIF ( particles(n)%x < pos_x_old  .AND. &
    312                    particles(n)%y > pos_y_old )  THEN
    313 
    314              t_index = 1
    315 
    316              DO  i = i1, i2, -1
    317                 xline      = i * dx - 0.5 * dx
    318                 t(t_index) = ( pos_x_old - xline ) / &
    319                              ( pos_x_old - particles(n)%x )
    320                 t_index    = t_index + 1
    321              ENDDO
    322 
    323              DO  j = j1, j2
    324                 yline      = j * dy + 0.5 * dy
    325                 t(t_index) = ( yline - pos_y_old ) / &
    326                              ( particles(n)%y - pos_y_old )
    327                 t_index    = t_index + 1
    328              ENDDO
    329 
    330              IF ( particles(n)%z < pos_z_old )  THEN
    331                 DO  k = k1, k2 , -1
    332                    zline      = k * dz
    333                    t(t_index) = ( pos_z_old - zline ) / &
    334                                 ( pos_z_old - particles(n)%z )
    335                    t_index    = t_index + 1
    336                 ENDDO
    337              ENDIF
    338              t_index_number = t_index - 1
    339 
    340 !
    341 !--          Sorting t
    342              inc = 1
    343              jr  = 1
    344 
    345              DO WHILE ( inc <= t_index_number )
    346                 inc = 3 * inc + 1
    347              ENDDO
    348 
    349              DO WHILE ( inc > 1 )
    350                 inc = inc / 3
    351                 DO  ir = inc+1, t_index_number
    352                    tmp_t = t(ir)
    353                    jr    = ir
    354                    DO WHILE ( t(jr-inc) > tmp_t )
    355                       t(jr) = t(jr-inc)
    356                       jr    = jr - inc
    357                       IF ( jr <= inc )