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

Last change on this file since 2716 was 2716, checked in by kanani, 4 years ago

Correction of "Former revisions" section

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