Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |