Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |