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

Last change on this file since 1929 was 1929, checked in by suehring, 8 years ago

several bugfixes in particle model and serial mode

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