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

Last change on this file since 2053 was 2001, checked in by knoop, 8 years ago

last commit documented

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