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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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