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

Last change on this file since 3046 was 2801, checked in by thiele, 3 years ago

Introduce particle transfer in nested models

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