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

Last change on this file since 3226 was 3189, checked in by Giersch, 6 years ago

Interpolation to particle position revised, indent revised, bugfix in calculation of the x/y indices for current particle position

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