Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |