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

Last change on this file since 3869 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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