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

Last change on this file since 2698 was 2698, checked in by suehring, 4 years ago

Particle reflections at downward-facing walls; revision of particle speed interpolations at walls; bugfixes in get_topography_index and in date constants

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