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

Last change on this file since 3436 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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