Home

Context Navigation

source: palm/tags/release-5.0/SOURCE/lpm_boundary_conds.f90 @ 4106

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |