Home

Context Navigation

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

Last change on this file since 2609 was 2606, checked in by schwenkel, 7 years ago
Modified particle box location and further changes in particle model
Property svn:keywords set to `Id`
File size: 25.5 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 2606 2017-11-10 10:36:31Z schwenkel $
27	! Changed particle box locations: center of particle box now coincides
28	! with scalar grid point of same index.
29	! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
30	! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
31	! lpm_sort -> lpm_sort_timeloop_done
32	!
33	! 2318 2017-07-20 17:27:44Z suehring
34	! Get topography top index via Function call
35	!
36	! 2317 2017-07-20 17:27:19Z suehring
37	!
38	! 2232 2017-05-30 17:47:52Z suehring
39	! Adjustments to new topography and surface concept
40	! Rename character range into location, as range is an intrinsic.
41	!
42	! 2000 2016-08-20 18:09:15Z knoop
43	! Forced header and separation lines into 80 columns
44	!
45	! 1929 2016-06-09 16:25:25Z suehring
46	! Rewritten wall reflection
47	!
48	! 1822 2016-04-07 07:49:42Z hoffmann
49	! Tails removed. Unused variables removed.
50	!
51	! 1682 2015-10-07 23:56:08Z knoop
52	! Code annotations made doxygen readable
53	!
54	! 1359 2014-04-11 17:15:14Z hoffmann
55	! New particle structure integrated.
56	! Kind definition added to all floating point numbers.
57	!
58	! 1320 2014-03-20 08:40:49Z raasch
59	! ONLY-attribute added to USE-statements,
60	! kind-parameters added to all INTEGER and REAL declaration statements,
61	! kinds are defined in new module kinds,
62	! revision history before 2012 removed,
63	! comment fields (!:) to be used for variable explanations added to
64	! all variable declaration statements
65	!
66	! 1036 2012-10-22 13:43:42Z raasch
67	! code put under GPL (PALM 3.9)
68	!
69	! 849 2012-03-15 10:35:09Z raasch
70	! routine renamed lpm_boundary_conds, bottom and top boundary conditions
71	! included (former part of advec_particles)
72	!
73	! 824 2012-02-17 09:09:57Z raasch
74	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3
75	!
76	! Initial version (2007/03/09)
77	!
78	! Description:
79	! ------------
80	!> Boundary conditions for the Lagrangian particles.
81	!> The routine consists of two different parts. One handles the bottom (flat)
82	!> and top boundary. In this part, also particles which exceeded their lifetime
83	!> are deleted.
84	!> The other part handles the reflection of particles from vertical walls.
85	!> This part was developed by Jin Zhang during 2006-2007.
86	!>
87	!> To do: Code structure for finding the t_index values and for checking the
88	!> ----- reflection conditions is basically the same for all four cases, so it
89	!> should be possible to further simplify/shorten it.
90	!>
91	!> THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
92	!> (see offset_ocean_*)
93	!------------------------------------------------------------------------------!
94	SUBROUTINE lpm_boundary_conds( location )
95
96
97	USE arrays_3d, &
98	ONLY: zu, zw
99
100	USE control_parameters, &
101	ONLY: dz, message_string, particle_maximum_age
102
103	USE cpulog, &
104	ONLY: cpu_log, log_point_s
105
106	USE grid_variables, &
107	ONLY: ddx, dx, ddy, dy
108
109	USE indices, &
110	ONLY: nxl, nxr, nyn, nys, nz, nzb
111
112	USE kinds
113
114	USE particle_attributes, &
115	ONLY: deleted_particles, ibc_par_b, ibc_par_t, number_of_particles, &
116	particles, particle_type, offset_ocean_nzt_m1, &
117	use_sgs_for_particles
118
119	USE pegrid
120
121	USE surface_mod, &
122	ONLY: get_topography_top_index
123
124	IMPLICIT NONE
125
126	CHARACTER (LEN=*) :: location !<
127
128	INTEGER(iwp) :: inc !< dummy for sorting algorithmus
129	INTEGER(iwp) :: ir !< dummy for sorting algorithmus
130	INTEGER(iwp) :: i1 !< grid index (x) of old particle position
131	INTEGER(iwp) :: i2 !< grid index (x) of current particle position
132	INTEGER(iwp) :: i3 !< grid index (x) of intermediate particle position
133	INTEGER(iwp) :: jr !< dummy for sorting algorithmus
134	INTEGER(iwp) :: j1 !< grid index (y) of old particle position
135	INTEGER(iwp) :: j2 !< grid index (x) of current particle position
136	INTEGER(iwp) :: j3 !< grid index (x) of intermediate particle position
137	INTEGER(iwp) :: k_wall !< vertical index of topography top
138	INTEGER(iwp) :: n !< particle number
139	INTEGER(iwp) :: t_index !< running index for intermediate particle timesteps in reflection algorithmus
140	INTEGER(iwp) :: t_index_number !< number of intermediate particle timesteps in reflection algorithmus
141	INTEGER(iwp) :: tmp_x !< dummy for sorting algorithmus
142	INTEGER(iwp) :: tmp_y !< dummy for sorting algorithmus
143
144	INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
145	INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions
146
147	LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary
148	LOGICAL :: cross_wall_y !< flag to check if particle reflection along y is necessary
149	LOGICAL :: downwards !< flag to check if particle reflection along z is necessary (only if particle move downwards)
150	LOGICAL :: reflect_x !< flag to check if particle is already reflected along x
151	LOGICAL :: reflect_y !< flag to check if particle is already reflected along y
152	LOGICAL :: reflect_z !< flag to check if particle is already reflected along z
153	LOGICAL :: tmp_reach_x !< dummy for sorting algorithmus
154	LOGICAL :: tmp_reach_y !< dummy for sorting algorithmus
155	LOGICAL :: tmp_reach_z !< dummy for sorting algorithmus
156	LOGICAL :: x_wall_reached !< flag to check if particle has already reached wall
157	LOGICAL :: y_wall_reached !< flag to check if particle has already reached wall
158
159	LOGICAL, DIMENSION(0:10) :: reach_x !< flag to check if particle is at a yz-wall
160	LOGICAL, DIMENSION(0:10) :: reach_y !< flag to check if particle is at a xz-wall
161	LOGICAL, DIMENSION(0:10) :: reach_z !< flag to check if particle is at a xy-wall
162
163	REAL(wp) :: dt_particle !< particle timestep
164	REAL(wp) :: dum !< dummy argument
165	REAL(wp) :: eps = 1E-10_wp !< security number to check if particle has reached a wall
166	REAL(wp) :: pos_x !< intermediate particle position (x)
167	REAL(wp) :: pos_x_old !< particle position (x) at previous particle timestep
168	REAL(wp) :: pos_y !< intermediate particle position (y)
169	REAL(wp) :: pos_y_old !< particle position (y) at previous particle timestep
170	REAL(wp) :: pos_z !< intermediate particle position (z)
171	REAL(wp) :: pos_z_old !< particle position (z) at previous particle timestep
172	REAL(wp) :: prt_x !< current particle position (x)
173	REAL(wp) :: prt_y !< current particle position (y)
174	REAL(wp) :: prt_z !< current particle position (z)
175	REAL(wp) :: t_old !< previous reflection time
176	REAL(wp) :: tmp_t !< dummy for sorting algorithmus
177	REAL(wp) :: xwall !< location of wall in x
178	REAL(wp) :: ywall !< location of wall in y
179	REAL(wp) :: zwall1 !< location of wall in z (old grid box)
180	REAL(wp) :: zwall2 !< location of wall in z (current grid box)
181	REAL(wp) :: zwall3 !< location of wall in z (old y, current x)
182	REAL(wp) :: zwall4 !< location of wall in z (current y, old x)
183
184	REAL(wp), DIMENSION(0:10) :: t !< reflection time
185
186
187	IF ( location == 'bottom/top' ) THEN
188
189	!
190	!-- Apply boundary conditions to those particles that have crossed the top or
191	!-- bottom boundary and delete those particles, which are older than allowed
192	DO n = 1, number_of_particles
193
194	!
195	!-- Stop if particles have moved further than the length of one
196	!-- PE subdomain (newly released particles have age = age_m!)
197	IF ( particles(n)%age /= particles(n)%age_m ) THEN
198	IF ( ABS(particles(n)%speed_x) > &
199	((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. &
200	ABS(particles(n)%speed_y) > &
201	((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN
202
203	WRITE( message_string, * ) 'particle too fast. n = ', n
204	CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
205	ENDIF
206	ENDIF
207
208	IF ( particles(n)%age > particle_maximum_age .AND. &
209	particles(n)%particle_mask ) &
210	THEN
211	particles(n)%particle_mask = .FALSE.
212	deleted_particles = deleted_particles + 1
213	ENDIF
214
215	IF ( particles(n)%z >= zu(nz) .AND. particles(n)%particle_mask ) THEN
216	IF ( ibc_par_t == 1 ) THEN
217	!
218	!-- Particle absorption
219	particles(n)%particle_mask = .FALSE.
220	deleted_particles = deleted_particles + 1
221	ELSEIF ( ibc_par_t == 2 ) THEN
222	!
223	!-- Particle reflection
224	particles(n)%z = 2.0_wp * zu(nz) - particles(n)%z
225	particles(n)%speed_z = -particles(n)%speed_z
226	IF ( use_sgs_for_particles .AND. &
227	particles(n)%rvar3 > 0.0_wp ) THEN
228	particles(n)%rvar3 = -particles(n)%rvar3
229	ENDIF
230	ENDIF
231	ENDIF
232
233	IF ( particles(n)%z < zw(0) .AND. particles(n)%particle_mask ) THEN
234	IF ( ibc_par_b == 1 ) THEN
235	!
236	!-- Particle absorption
237	particles(n)%particle_mask = .FALSE.
238	deleted_particles = deleted_particles + 1
239	ELSEIF ( ibc_par_b == 2 ) THEN
240	!
241	!-- Particle reflection
242	particles(n)%z = 2.0_wp * zw(0) - particles(n)%z
243	particles(n)%speed_z = -particles(n)%speed_z
244	IF ( use_sgs_for_particles .AND. &
245	particles(n)%rvar3 < 0.0_wp ) THEN
246	particles(n)%rvar3 = -particles(n)%rvar3
247	ENDIF
248	ENDIF
249	ENDIF
250	ENDDO
251
252	ELSEIF ( location == 'walls' ) THEN
253
254
255	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
256
257	DO n = 1, number_of_particles
258	!
259	!-- Recalculate particle timestep
260	dt_particle = particles(n)%age - particles(n)%age_m
261	!
262	!-- Obtain x/y indices for current particle position
263	i2 = particles(n)%x * ddx
264	j2 = particles(n)%y * ddy
265	!
266	!-- Save current particle positions
267	prt_x = particles(n)%x
268	prt_y = particles(n)%y
269	prt_z = particles(n)%z
270	!
271	!-- Recalculate old particle positions
272	pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
273	pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
274	pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
275	!
276	!-- Obtain x/y indices for old particle positions
277	i1 = pos_x_old * ddx
278	j1 = pos_y_old * ddy
279	!
280	!-- Determine horizontal as well as vertical walls at which particle can
281	!-- be potentially reflected.
282	!-- Start with walls aligned in yz layer.
283	!-- Wall to the right
284	IF ( prt_x > pos_x_old ) THEN
285	xwall = ( i1 + 0.5_wp ) * dx
286	!
287	!-- Wall to the left
288	ELSE
289	xwall = ( i1 - 0.5_wp ) * dx
290	ENDIF
291	!
292	!-- Walls aligned in xz layer
293	!-- Wall to the north
294	IF ( prt_y > pos_y_old ) THEN
295	ywall = ( j1 + 0.5_wp ) * dy
296	!-- Wall to the south
297	ELSE
298	ywall = ( j1 - 0.5_wp ) * dy
299	ENDIF
300	!
301	!-- Walls aligned in xy layer at which particle can be possiblly reflected.
302	!-- The construct of MERGE and BTEST is used to determine the topography-
303	!-- top index (former nzb_s_inner).
304	zwall1 = zw( get_topography_top_index( j2, i2, 's' ) )
305	zwall2 = zw( get_topography_top_index( j1, i1, 's' ) )
306	zwall3 = zw( get_topography_top_index( j1, i2, 's' ) )
307	zwall4 = zw( get_topography_top_index( j2, i1, 's' ) )
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 = get_topography_top_index( j3, i3, 's' )
493	IF ( ABS( pos_x - xwall ) < eps .AND. &
494	pos_z <= zw(k_wall) .AND. &
495	reach_x(t_index) .AND. &
496	.NOT. reflect_x ) THEN
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 do not exactly matches 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	!-- Change also 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 particle do not crosses 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
525	!
526	!-- Check if a particle needs to be reflected at any xz-wall. If
527	!-- necessary, carry out reflection. At first, determine index of
528	!-- topography top at (j3,i3)
529	k_wall = get_topography_top_index( j3, i3, 's' )
530	IF ( ABS( pos_y - ywall ) < eps .AND. &
531	pos_z <= zw(k_wall) .AND. &
532	reach_y(t_index) .AND. &
533	.NOT. reflect_y ) THEN
534
535	pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ), &
536	MAX( 2.0_wp * ywall - pos_y, ywall ), &
537	particles(n)%y > ywall )
538
539	particles(n)%speed_y = - particles(n)%speed_y
540	particles(n)%rvar2 = - particles(n)%rvar2
541
542	reflect_y = .TRUE.
543	y_ind(t_index:t_index_number) = j1
544
545	ELSEIF ( y_wall_reached .AND. .NOT. reflect_y ) THEN
546	y_ind(t_index:t_index_number) = j2
547	ENDIF
548	!
549	!-- Check if a particle needs to be reflected at any xy-wall. If
550	!-- necessary, carry out reflection.
551	IF ( downwards .AND. reach_z(t_index) .AND. &
552	.NOT. reflect_z ) THEN
553	!
554	!-- Determine index of topography top at (j3,i3) and chick if
555	!-- particle is below.
556	k_wall = get_topography_top_index( j3, i3, 's' )
557	IF ( pos_z - zw(k_wall) < eps ) THEN
558
559	pos_z = MAX( 2.0_wp * zw(k_wall) - pos_z, &
560	zw(k_wall) )
561
562	particles(n)%speed_z = - particles(n)%speed_z
563	particles(n)%rvar3 = - particles(n)%rvar3
564
565	reflect_z = .TRUE.
566
567	ENDIF
568
569	ENDIF
570	!
571	!-- Swap time
572	t_old = t(t_index)
573
574	ENDDO
575	!
576	!-- If a particle was reflected, calculate final position from last
577	!-- intermediate position.
578	IF ( reflect_x .OR. reflect_y .OR. reflect_z ) THEN
579
580	particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle &
581	* particles(n)%speed_x
582	particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle &
583	* particles(n)%speed_y
584	particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle &
585	* particles(n)%speed_z
586
587	ENDIF
588
589	ENDIF
590
591	ENDDO
592
593	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
594
595	ENDIF
596
597	END SUBROUTINE lpm_boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |