Home

Context Navigation

advec_particles.f90 @ 667

Last change on this file since 667 was 667, checked in by suehring, 14 years ago

summary:

Gryschka:

Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
Exchange of u and v from ocean surface to atmosphere surface
Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
Inflow turbulence is now defined by flucuations around spanwise mean
Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

New advection added ( Wicker and Skamarock 5th order ), therefore:
- New module advec_ws.f90
- Modified exchange of ghost boundaries.
- Modified evaluation of turbulent fluxes
- New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90

Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90

Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90

Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90

When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90

Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90

Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90

Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90

Steering parameter dissipation_control added in inipar.

Makefile

Module advec_ws added.

Modules

Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90

New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90

Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90

Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90

Output of total array size was adapted to nbgp.

data_output_2d.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90

Calls of exchange_horiz are modified.

diffusion_e.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90

Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90

Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90

determination of target_id's moved to init_pegrid

init_pt_anomaly.f90

Call of exchange_horiz are modified.

init_rankine.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90

Output of advection scheme.

poismg.f90

Calls of exchange_horiz are modified.

prandtl_fluxes.f90

Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90

+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90

Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90

Allocation with nbgp.

wall_fluxes.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90

Array bounds and nx, ny adapted with nbgp.

sor.f90

bug removed in declaration of ddzw(), nz replaced by nzt+1

Property svn:keywords set to Id

File size: 171.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/SOURCE/advec_particles.f90 @ 667

Download in other formats: