Home

Context Navigation

source: palm/trunk/SOURCE/lpm_boundary_conds.f90 @ 2284

Last change on this file since 2284 was 2233, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 27.0 KB

Line
1	!> @file lpm_boundary_conds.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: lpm_boundary_conds.f90 2233 2017-05-30 18:08:54Z witha $
27	!
28	! 2232 2017-05-30 17:47:52Z suehring
29	! Adjustments to new topography and surface concept
30	! Rename character range into location, as range is an intrinsic.
31	!
32	! 2000 2016-08-20 18:09:15Z knoop
33	! Forced header and separation lines into 80 columns
34	!
35	! 1929 2016-06-09 16:25:25Z suehring
36	! Rewritten wall reflection
37	!
38	! 1822 2016-04-07 07:49:42Z hoffmann
39	! Tails removed. Unused variables removed.
40	!
41	! 1682 2015-10-07 23:56:08Z knoop
42	! Code annotations made doxygen readable
43	!
44	! 1359 2014-04-11 17:15:14Z hoffmann
45	! New particle structure integrated.
46	! Kind definition added to all floating point numbers.
47	!
48	! 1320 2014-03-20 08:40:49Z raasch
49	! ONLY-attribute added to USE-statements,
50	! kind-parameters added to all INTEGER and REAL declaration statements,
51	! kinds are defined in new module kinds,
52	! revision history before 2012 removed,
53	! comment fields (!:) to be used for variable explanations added to
54	! all variable declaration statements
55	!
56	! 1036 2012-10-22 13:43:42Z raasch
57	! code put under GPL (PALM 3.9)
58	!
59	! 849 2012-03-15 10:35:09Z raasch
60	! routine renamed lpm_boundary_conds, bottom and top boundary conditions
61	! included (former part of advec_particles)
62	!
63	! 824 2012-02-17 09:09:57Z raasch
64	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3
65	!
66	! Initial version (2007/03/09)
67	!
68	! Description:
69	! ------------
70	!> Boundary conditions for the Lagrangian particles.
71	!> The routine consists of two different parts. One handles the bottom (flat)
72	!> and top boundary. In this part, also particles which exceeded their lifetime
73	!> are deleted.
74	!> The other part handles the reflection of particles from vertical walls.
75	!> This part was developed by Jin Zhang during 2006-2007.
76	!>
77	!> To do: Code structure for finding the t_index values and for checking the
78	!> ----- reflection conditions is basically the same for all four cases, so it
79	!> should be possible to further simplify/shorten it.
80	!>
81	!> THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
82	!> (see offset_ocean_*)
83	!------------------------------------------------------------------------------!
84	SUBROUTINE lpm_boundary_conds( location )
85
86
87	USE arrays_3d, &
88	ONLY: zu, zw
89
90	USE control_parameters, &
91	ONLY: dz, message_string, particle_maximum_age
92
93	USE cpulog, &
94	ONLY: cpu_log, log_point_s
95
96	USE grid_variables, &
97	ONLY: ddx, dx, ddy, dy
98
99	USE indices, &
100	ONLY: nxl, nxr, nyn, nys, nz, nzb, nzb_max, wall_flags_0
101
102	USE kinds
103
104	USE particle_attributes, &
105	ONLY: deleted_particles, ibc_par_b, ibc_par_t, number_of_particles, &
106	particles, particle_type, offset_ocean_nzt_m1, &
107	use_sgs_for_particles
108
109	USE pegrid
110
111	IMPLICIT NONE
112
113	CHARACTER (LEN=*) :: location !<
114
115	INTEGER(iwp) :: inc !< dummy for sorting algorithmus
116	INTEGER(iwp) :: ir !< dummy for sorting algorithmus
117	INTEGER(iwp) :: i1 !< grid index (x) of old particle position
118	INTEGER(iwp) :: i2 !< grid index (x) of current particle position
119	INTEGER(iwp) :: i3 !< grid index (x) of intermediate particle position
120	INTEGER(iwp) :: jr !< dummy for sorting algorithmus
121	INTEGER(iwp) :: j1 !< grid index (y) of old particle position
122	INTEGER(iwp) :: j2 !< grid index (x) of current particle position
123	INTEGER(iwp) :: j3 !< grid index (x) of intermediate particle position
124	INTEGER(iwp) :: k_wall !< vertical index of topography top
125	INTEGER(iwp) :: n !< particle number
126	INTEGER(iwp) :: t_index !< running index for intermediate particle timesteps in reflection algorithmus
127	INTEGER(iwp) :: t_index_number !< number of intermediate particle timesteps in reflection algorithmus
128	INTEGER(iwp) :: tmp_x !< dummy for sorting algorithmus
129	INTEGER(iwp) :: tmp_y !< dummy for sorting algorithmus
130
131	INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
132	INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions
133
134	LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary
135	LOGICAL :: cross_wall_y !< flag to check if particle reflection along y is necessary
136	LOGICAL :: downwards !< flag to check if particle reflection along z is necessary (only if particle move downwards)
137	LOGICAL :: reflect_x !< flag to check if particle is already reflected along x
138	LOGICAL :: reflect_y !< flag to check if particle is already reflected along y
139	LOGICAL :: reflect_z !< flag to check if particle is already reflected along z
140	LOGICAL :: tmp_reach_x !< dummy for sorting algorithmus
141	LOGICAL :: tmp_reach_y !< dummy for sorting algorithmus
142	LOGICAL :: tmp_reach_z !< dummy for sorting algorithmus
143	LOGICAL :: x_wall_reached !< flag to check if particle has already reached wall
144	LOGICAL :: y_wall_reached !< flag to check if particle has already reached wall
145
146	LOGICAL, DIMENSION(0:10) :: reach_x !< flag to check if particle is at a yz-wall
147	LOGICAL, DIMENSION(0:10) :: reach_y !< flag to check if particle is at a xz-wall
148	LOGICAL, DIMENSION(0:10) :: reach_z !< flag to check if particle is at a xy-wall
149
150	REAL(wp) :: dt_particle !< particle timestep
151	REAL(wp) :: dum !< dummy argument
152	REAL(wp) :: eps = 1E-10_wp !< security number to check if particle has reached a wall
153	REAL(wp) :: pos_x !< intermediate particle position (x)
154	REAL(wp) :: pos_x_old !< particle position (x) at previous particle timestep
155	REAL(wp) :: pos_y !< intermediate particle position (y)
156	REAL(wp) :: pos_y_old !< particle position (y) at previous particle timestep
157	REAL(wp) :: pos_z !< intermediate particle position (z)
158	REAL(wp) :: pos_z_old !< particle position (z) at previous particle timestep
159	REAL(wp) :: prt_x !< current particle position (x)
160	REAL(wp) :: prt_y !< current particle position (y)
161	REAL(wp) :: prt_z !< current particle position (z)
162	REAL(wp) :: t_old !< previous reflection time
163	REAL(wp) :: tmp_t !< dummy for sorting algorithmus
164	REAL(wp) :: xwall !< location of wall in x
165	REAL(wp) :: ywall !< location of wall in y
166	REAL(wp) :: zwall1 !< location of wall in z (old grid box)
167	REAL(wp) :: zwall2 !< location of wall in z (current grid box)
168	REAL(wp) :: zwall3 !< location of wall in z (old y, current x)
169	REAL(wp) :: zwall4 !< location of wall in z (current y, old x)
170
171	REAL(wp), DIMENSION(0:10) :: t !< reflection time
172
173
174	IF ( location == 'bottom/top' ) THEN
175
176	!
177	!-- Apply boundary conditions to those particles that have crossed the top or
178	!-- bottom boundary and delete those particles, which are older than allowed
179	DO n = 1, number_of_particles
180
181	!
182	!-- Stop if particles have moved further than the length of one
183	!-- PE subdomain (newly released particles have age = age_m!)
184	IF ( particles(n)%age /= particles(n)%age_m ) THEN
185	IF ( ABS(particles(n)%speed_x) > &
186	((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. &
187	ABS(particles(n)%speed_y) > &
188	((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN
189
190	WRITE( message_string, * ) 'particle too fast. n = ', n
191	CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
192	ENDIF
193	ENDIF
194
195	IF ( particles(n)%age > particle_maximum_age .AND. &
196	particles(n)%particle_mask ) &
197	THEN
198	particles(n)%particle_mask = .FALSE.
199	deleted_particles = deleted_particles + 1
200	ENDIF
201
202	IF ( particles(n)%z >= zu(nz) .AND. particles(n)%particle_mask ) THEN
203	IF ( ibc_par_t == 1 ) THEN
204	!
205	!-- Particle absorption
206	particles(n)%particle_mask = .FALSE.
207	deleted_particles = deleted_particles + 1
208	ELSEIF ( ibc_par_t == 2 ) THEN
209	!
210	!-- Particle reflection
211	particles(n)%z = 2.0_wp * zu(nz) - particles(n)%z
212	particles(n)%speed_z = -particles(n)%speed_z
213	IF ( use_sgs_for_particles .AND. &
214	particles(n)%rvar3 > 0.0_wp ) THEN
215	particles(n)%rvar3 = -particles(n)%rvar3
216	ENDIF
217	ENDIF
218	ENDIF
219
220	IF ( particles(n)%z < zw(0) .AND. particles(n)%particle_mask ) THEN
221	IF ( ibc_par_b == 1 ) THEN
222	!
223	!-- Particle absorption
224	particles(n)%particle_mask = .FALSE.
225	deleted_particles = deleted_particles + 1
226	ELSEIF ( ibc_par_b == 2 ) THEN
227	!
228	!-- Particle reflection
229	particles(n)%z = 2.0_wp * zw(0) - particles(n)%z
230	particles(n)%speed_z = -particles(n)%speed_z
231	IF ( use_sgs_for_particles .AND. &
232	particles(n)%rvar3 < 0.0_wp ) THEN
233	particles(n)%rvar3 = -particles(n)%rvar3
234	ENDIF
235	ENDIF
236	ENDIF
237	ENDDO
238
239	ELSEIF ( location == 'walls' ) THEN
240
241
242	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
243
244	DO n = 1, number_of_particles
245	!
246	!-- Recalculate particle timestep
247	dt_particle = particles(n)%age - particles(n)%age_m
248	!
249	!-- Obtain x/y indices for current particle position
250	i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
251	j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
252	!
253	!-- Save current particle positions
254	prt_x = particles(n)%x
255	prt_y = particles(n)%y
256	prt_z = particles(n)%z
257	!
258	!-- Recalculate old particle positions
259	pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
260	pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
261	pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
262	!
263	!-- Obtain x/y indices for old particle positions
264	i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
265	j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
266	!
267	!-- Determine horizontal as well as vertical walls at which particle can
268	!-- be potentially reflected.
269	!-- Start with walls aligned in yz layer.
270	!-- Wall to the right
271	IF ( prt_x > pos_x_old ) THEN
272	xwall = ( i1 + 0.5_wp ) * dx
273	!
274	!-- Wall to the left
275	ELSE
276	xwall = ( i1 - 0.5_wp ) * dx
277	ENDIF
278	!
279	!-- Walls aligned in xz layer
280	!-- Wall to the north
281	IF ( prt_y > pos_y_old ) THEN
282	ywall = ( j1 + 0.5_wp ) * dy
283	!-- Wall to the south
284	ELSE
285	ywall = ( j1 - 0.5_wp ) * dy
286	ENDIF
287	!
288	!-- Walls aligned in xy layer at which particle can be possiblly reflected.
289	!-- The construct of MERGE and BTEST is used to determine the topography-
290	!-- top index (former nzb_s_inner).
291	zwall1 = zw( MAXLOC( &
292	MERGE( 1, 0, &
293	BTEST( wall_flags_0(nzb:nzb_max,j2,i2), 12 ) &
294	), DIM = 1 &
295	) - 1 )
296	zwall2 = zw( MAXLOC( &
297	MERGE( 1, 0, &
298	BTEST( wall_flags_0(nzb:nzb_max,j1,i1), 12 ) &
299	), DIM = 1 &
300	) - 1 )
301	zwall3 = zw( MAXLOC( &
302	MERGE( 1, 0, &
303	BTEST( wall_flags_0(nzb:nzb_max,j1,i2), 12 ) &
304	), DIM = 1 &
305	) - 1 )
306	zwall4 = zw( MAXLOC( &
307	MERGE( 1, 0, &
308	BTEST( wall_flags_0(nzb:nzb_max,j2,i1), 12 ) &
309	), DIM = 1 &
310	) - 1 )
311	!
312	!-- Initialize flags to check if particle reflection is necessary
313	downwards = .FALSE.
314	cross_wall_x = .FALSE.
315	cross_wall_y = .FALSE.
316	!
317	!-- Initialize flags to check if a wall is reached
318	reach_x = .FALSE.
319	reach_y = .FALSE.
320	reach_z = .FALSE.
321	!
322	!-- Initialize flags to check if a particle was already reflected
323	reflect_x = .FALSE.
324	reflect_y = .FALSE.
325	reflect_z = .FALSE.
326	!
327	!-- Initialize flags to check if a vertical wall is already crossed.
328	!-- ( Required to obtain correct indices. )
329	x_wall_reached = .FALSE.
330	y_wall_reached = .FALSE.
331	!
332	!-- Initialize time array
333	t = 0.0_wp
334	!
335	!-- Check if particle can reach any wall. This case, calculate the
336	!-- fractional time needed to reach this wall. Store this fractional
337	!-- timestep in array t. Moreover, store indices for these grid
338	!-- boxes where the respective wall belongs to.
339	!-- Start with x-direction.
340	t_index = 1
341	t(t_index) = ( xwall - pos_x_old ) &
342	/ MERGE( MAX( prt_x - pos_x_old, 1E-30_wp ), &
343	MIN( prt_x - pos_x_old, -1E-30_wp ), &
344	prt_x > pos_x_old )
345	x_ind(t_index) = i2
346	y_ind(t_index) = j1
347	reach_x(t_index) = .TRUE.
348	reach_y(t_index) = .FALSE.
349	reach_z(t_index) = .FALSE.
350	!
351	!-- Store these values only if particle really reaches any wall. t must
352	!-- be in a interval between [0:1].
353	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN
354	t_index = t_index + 1
355	cross_wall_x = .TRUE.
356	ENDIF
357	!
358	!-- y-direction
359	t(t_index) = ( ywall - pos_y_old ) &
360	/ MERGE( MAX( prt_y - pos_y_old, 1E-30_wp ), &
361	MIN( prt_y - pos_y_old, -1E-30_wp ), &
362	prt_y > pos_y_old )
363	x_ind(t_index) = i1
364	y_ind(t_index) = j2
365	reach_x(t_index) = .FALSE.
366	reach_y(t_index) = .TRUE.
367	reach_z(t_index) = .FALSE.
368	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN
369	t_index = t_index + 1
370	cross_wall_y = .TRUE.
371	ENDIF
372	!
373	!-- z-direction
374	!-- At first, check if particle moves downwards. Only in this case a
375	!-- particle can be reflected vertically.
376	IF ( prt_z < pos_z_old ) THEN
377
378	downwards = .TRUE.
379	dum = 1.0_wp / MERGE( MAX( prt_z - pos_z_old, 1E-30_wp ), &
380	MIN( prt_z - pos_z_old, -1E-30_wp ), &
381	prt_z > pos_z_old )
382
383	t(t_index) = ( zwall1 - pos_z_old ) * dum
384	x_ind(t_index) = i2
385	y_ind(t_index) = j2
386	reach_x(t_index) = .FALSE.
387	reach_y(t_index) = .FALSE.
388	reach_z(t_index) = .TRUE.
389	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) &
390	t_index = t_index + 1
391
392	reach_x(t_index) = .FALSE.
393	reach_y(t_index) = .FALSE.
394	reach_z(t_index) = .TRUE.
395	t(t_index) = ( zwall2 - pos_z_old ) * dum
396	x_ind(t_index) = i1
397	y_ind(t_index) = j1
398	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) &
399	t_index = t_index + 1
400
401	reach_x(t_index) = .FALSE.
402	reach_y(t_index) = .FALSE.
403	reach_z(t_index) = .TRUE.
404	t(t_index) = ( zwall3 - pos_z_old ) * dum
405	x_ind(t_index) = i2
406	y_ind(t_index) = j1
407	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) &
408	t_index = t_index + 1
409
410	reach_x(t_index) = .FALSE.
411	reach_y(t_index) = .FALSE.
412	reach_z(t_index) = .TRUE.
413	t(t_index) = ( zwall4 - pos_z_old ) * dum
414	x_ind(t_index) = i1
415	y_ind(t_index) = j2
416	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) &
417	t_index = t_index + 1
418
419	END IF
420	t_index_number = t_index - 1
421	!
422	!-- Carry out reflection only if particle reaches any wall
423	IF ( cross_wall_x .OR. cross_wall_y .OR. downwards ) THEN
424	!
425	!-- Sort fractional timesteps in ascending order. Also sort the
426	!-- corresponding indices and flag according to the time interval a
427	!-- particle reaches the respective wall.
428	inc = 1
429	jr = 1
430	DO WHILE ( inc <= t_index_number )
431	inc = 3 * inc + 1
432	ENDDO
433
434	DO WHILE ( inc > 1 )
435	inc = inc / 3
436	DO ir = inc+1, t_index_number
437	tmp_t = t(ir)
438	tmp_x = x_ind(ir)
439	tmp_y = y_ind(ir)
440	tmp_reach_x = reach_x(ir)
441	tmp_reach_y = reach_y(ir)
442	tmp_reach_z = reach_z(ir)
443	jr = ir
444	DO WHILE ( t(jr-inc) > tmp_t )
445	t(jr) = t(jr-inc)
446	x_ind(jr) = x_ind(jr-inc)
447	y_ind(jr) = y_ind(jr-inc)
448	reach_x(jr) = reach_x(jr-inc)
449	reach_y(jr) = reach_y(jr-inc)
450	reach_z(jr) = reach_z(jr-inc)
451	jr = jr - inc
452	IF ( jr <= inc ) EXIT
453	ENDDO
454	t(jr) = tmp_t
455	x_ind(jr) = tmp_x
456	y_ind(jr) = tmp_y
457	reach_x(jr) = tmp_reach_x
458	reach_y(jr) = tmp_reach_y
459	reach_z(jr) = tmp_reach_z
460	ENDDO
461	ENDDO
462	!
463	!-- Initialize temporary particle positions
464	pos_x = pos_x_old
465	pos_y = pos_y_old
466	pos_z = pos_z_old
467	!
468	!-- Loop over all times a particle possibly moves into a new grid box
469	t_old = 0.0_wp
470	DO t_index = 1, t_index_number
471	!
472	!-- Calculate intermediate particle position according to the
473	!-- timesteps a particle reaches any wall.
474	pos_x = pos_x + ( t(t_index) - t_old ) * dt_particle &
475	* particles(n)%speed_x
476	pos_y = pos_y + ( t(t_index) - t_old ) * dt_particle &
477	* particles(n)%speed_y
478	pos_z = pos_z + ( t(t_index) - t_old ) * dt_particle &
479	* particles(n)%speed_z
480	!
481	!-- Obtain x/y grid indices for intermediate particle position from
482	!-- sorted index array
483	i3 = x_ind(t_index)
484	j3 = y_ind(t_index)
485	!
486	!-- Check which wall is already reached
487	IF ( .NOT. x_wall_reached ) x_wall_reached = reach_x(t_index)
488	IF ( .NOT. y_wall_reached ) y_wall_reached = reach_y(t_index)
489	!
490	!-- Check if a particle needs to be reflected at any yz-wall. If
491	!-- necessary, carry out reflection. Please note, a security
492	!-- constant is required, as the particle position do not
493	!-- necessarily exactly match the wall location due to rounding
494	!-- errors. At first, determine index of topography top at (j3,i3)
495	k_wall = MAXLOC( &
496	MERGE( 1, 0, &
497	BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 ) &
498	), DIM = 1 &
499	) - 1
500	IF ( ABS( pos_x - xwall ) < eps .AND. &
501	pos_z <= zw(k_wall) .AND. &
502	reach_x(t_index) .AND. &
503	.NOT. reflect_x ) THEN
504	!
505	!-- Reflection in x-direction.
506	!-- Ensure correct reflection by MIN/MAX functions, depending on
507	!-- direction of particle transport.
508	!-- Due to rounding errors pos_x do not exactly matches the wall
509	!-- location, leading to erroneous reflection.
510	pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ), &
511	MAX( 2.0_wp * xwall - pos_x, xwall ), &
512	particles(n)%x > xwall )
513	!
514	!-- Change sign of particle speed
515	particles(n)%speed_x = - particles(n)%speed_x
516	!
517	!-- Change also sign of subgrid-scale particle speed
518	particles(n)%rvar1 = - particles(n)%rvar1
519	!
520	!-- Set flag that reflection along x is already done
521	reflect_x = .TRUE.
522	!
523	!-- As particle do not crosses any further yz-wall during
524	!-- this timestep, set further x-indices to the current one.
525	x_ind(t_index:t_index_number) = i1
526	!
527	!-- If particle already reached the wall but was not reflected,
528	!-- set further x-indices to the new one.
529	ELSEIF ( x_wall_reached .AND. .NOT. reflect_x ) THEN
530	x_ind(t_index:t_index_number) = i2
531	ENDIF
532	!
533	!-- Check if a particle needs to be reflected at any xz-wall. If
534	!-- necessary, carry out reflection. At first, determine index of
535	!-- topography top at (j3,i3)
536	k_wall = MAXLOC( &
537	MERGE( 1, 0, &
538	BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 ) &
539	), DIM = 1 &
540	) - 1
541	IF ( ABS( pos_y - ywall ) < eps .AND. &
542	pos_z <= zw(k_wall) .AND. &
543	reach_y(t_index) .AND. &
544	.NOT. reflect_y ) THEN
545
546	pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ), &
547	MAX( 2.0_wp * ywall - pos_y, ywall ), &
548	particles(n)%y > ywall )
549
550	particles(n)%speed_y = - particles(n)%speed_y
551	particles(n)%rvar2 = - particles(n)%rvar2
552
553	reflect_y = .TRUE.
554	y_ind(t_index:t_index_number) = j1
555
556	ELSEIF ( y_wall_reached .AND. .NOT. reflect_y ) THEN
557	y_ind(t_index:t_index_number) = j2
558	ENDIF
559	!
560	!-- Check if a particle needs to be reflected at any xy-wall. If
561	!-- necessary, carry out reflection.
562	IF ( downwards .AND. reach_z(t_index) .AND. &
563	.NOT. reflect_z ) THEN
564	!
565	!-- Determine index of topography top at (j3,i3) and chick if
566	!-- particle is below.
567	k_wall = MAXLOC( &
568	MERGE( 1, 0, &
569	BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 )&
570	), DIM = 1 &
571	) - 1
572	IF ( pos_z - zw(k_wall) < eps ) THEN
573
574	pos_z = MAX( 2.0_wp * zw(k_wall) - pos_z, &
575	zw(k_wall) )
576
577	particles(n)%speed_z = - particles(n)%speed_z
578	particles(n)%rvar3 = - particles(n)%rvar3
579
580	reflect_z = .TRUE.
581
582	ENDIF
583
584	ENDIF
585	!
586	!-- Swap time
587	t_old = t(t_index)
588
589	ENDDO
590	!
591	!-- If a particle was reflected, calculate final position from last
592	!-- intermediate position.
593	IF ( reflect_x .OR. reflect_y .OR. reflect_z ) THEN
594
595	particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle &
596	* particles(n)%speed_x
597	particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle &
598	* particles(n)%speed_y
599	particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle &
600	* particles(n)%speed_z
601
602	ENDIF
603
604	ENDIF
605
606	ENDDO
607
608	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
609
610	ENDIF
611
612	END SUBROUTINE lpm_boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |