Home

Context Navigation

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

Last change on this file since 2232 was 2232, checked in by suehring, 7 years ago
Adjustments according new topography and surface-modelling concept implemented
Property svn:keywords set to `Id`
File size: 27.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |