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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

  • Property svn:keywords set to Id
File size: 27.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
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! Adjustments to new topography and surface concept
23! Rename character range into location, as range is an intrinsic.
24!
25! Former revisions:
26! -----------------
27! $Id: lpm_boundary_conds.f90 2232 2017-05-30 17:47:52Z suehring $
28!
29! 2000 2016-08-20 18:09:15Z knoop
30! Forced header and separation lines into 80 columns
31!
32! 1929 2016-06-09 16:25:25Z suehring
33! Rewritten wall reflection
34!
35! 1822 2016-04-07 07:49:42Z hoffmann
36! Tails removed. Unused variables removed.
37!
38! 1682 2015-10-07 23:56:08Z knoop
39! Code annotations made doxygen readable
40!
41! 1359 2014-04-11 17:15:14Z hoffmann
42! New particle structure integrated.
43! Kind definition added to all floating point numbers.
44!
45! 1320 2014-03-20 08:40:49Z raasch
46! ONLY-attribute added to USE-statements,
47! kind-parameters added to all INTEGER and REAL declaration statements,
48! kinds are defined in new module kinds,
49! revision history before 2012 removed,
50! comment fields (!:) to be used for variable explanations added to
51! all variable declaration statements
52!
53! 1036 2012-10-22 13:43:42Z raasch
54! code put under GPL (PALM 3.9)
55!
56! 849 2012-03-15 10:35:09Z raasch
57! routine renamed lpm_boundary_conds, bottom and top boundary conditions
58! included (former part of advec_particles)
59!
60! 824 2012-02-17 09:09:57Z raasch
61! particle attributes speed_x|y|z_sgs renamed rvar1|2|3
62!
63! Initial version (2007/03/09)
64!
65! Description:
66! ------------
67!> Boundary conditions for the Lagrangian particles.
68!> The routine consists of two different parts. One handles the bottom (flat)
69!> and top boundary. In this part, also particles which exceeded their lifetime
70!> are deleted.
71!> The other part handles the reflection of particles from vertical walls.
72!> This part was developed by Jin Zhang during 2006-2007.
73!>
74!> To do: Code structure for finding the t_index values and for checking the
75!> -----  reflection conditions is basically the same for all four cases, so it
76!>        should be possible to further simplify/shorten it.
77!>
78!> THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
79!> (see offset_ocean_*)
80!------------------------------------------------------------------------------!
81 SUBROUTINE lpm_boundary_conds( location )
82 
83
84    USE arrays_3d,                                                             &
85        ONLY:  zu, zw
86
87    USE control_parameters,                                                    &
88        ONLY:  dz, message_string, particle_maximum_age
89
90    USE cpulog,                                                                &
91        ONLY:  cpu_log, log_point_s
92
93    USE grid_variables,                                                        &
94        ONLY:  ddx, dx, ddy, dy
95
96    USE indices,                                                               &
97        ONLY:  nxl, nxr, nyn, nys, nz, nzb, nzb_max, wall_flags_0
98
99    USE kinds
100
101    USE particle_attributes,                                                   &
102        ONLY:  deleted_particles, ibc_par_b, ibc_par_t, number_of_particles,   &
103               particles, particle_type, offset_ocean_nzt_m1,                  &
104               use_sgs_for_particles
105
106    USE pegrid
107
108    IMPLICIT NONE
109
110    CHARACTER (LEN=*) ::  location     !<
111   
112    INTEGER(iwp) ::  inc            !< dummy for sorting algorithmus
113    INTEGER(iwp) ::  ir             !< dummy for sorting algorithmus
114    INTEGER(iwp) ::  i1             !< grid index (x) of old particle position
115    INTEGER(iwp) ::  i2             !< grid index (x) of current particle position
116    INTEGER(iwp) ::  i3             !< grid index (x) of intermediate particle position
117    INTEGER(iwp) ::  jr             !< dummy for sorting algorithmus
118    INTEGER(iwp) ::  j1             !< grid index (y) of old particle position
119    INTEGER(iwp) ::  j2             !< grid index (x) of current particle position
120    INTEGER(iwp) ::  j3             !< grid index (x) of intermediate particle position
121    INTEGER(iwp) ::  k_wall         !< vertical index of topography top
122    INTEGER(iwp) ::  n              !< particle number
123    INTEGER(iwp) ::  t_index        !< running index for intermediate particle timesteps in reflection algorithmus
124    INTEGER(iwp) ::  t_index_number !< number of intermediate particle timesteps in reflection algorithmus
125    INTEGER(iwp) ::  tmp_x          !< dummy for sorting algorithmus
126    INTEGER(iwp) ::  tmp_y          !< dummy for sorting algorithmus
127
128    INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
129    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions
130   
131    LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
132    LOGICAL  ::  cross_wall_y    !< flag to check if particle reflection along y is necessary
133    LOGICAL  ::  downwards       !< flag to check if particle reflection along z is necessary (only if particle move downwards)
134    LOGICAL  ::  reflect_x       !< flag to check if particle is already reflected along x
135    LOGICAL  ::  reflect_y       !< flag to check if particle is already reflected along y
136    LOGICAL  ::  reflect_z       !< flag to check if particle is already reflected along z
137    LOGICAL  ::  tmp_reach_x     !< dummy for sorting algorithmus
138    LOGICAL  ::  tmp_reach_y     !< dummy for sorting algorithmus
139    LOGICAL  ::  tmp_reach_z     !< dummy for sorting algorithmus
140    LOGICAL  ::  x_wall_reached  !< flag to check if particle has already reached wall
141    LOGICAL  ::  y_wall_reached  !< flag to check if particle has already reached wall
142
143    LOGICAL, DIMENSION(0:10) ::  reach_x  !< flag to check if particle is at a yz-wall
144    LOGICAL, DIMENSION(0:10) ::  reach_y  !< flag to check if particle is at a xz-wall
145    LOGICAL, DIMENSION(0:10) ::  reach_z  !< flag to check if particle is at a xy-wall
146
147    REAL(wp) ::  dt_particle    !< particle timestep
148    REAL(wp) ::  dum            !< dummy argument
149    REAL(wp) ::  eps = 1E-10_wp !< security number to check if particle has reached a wall
150    REAL(wp) ::  pos_x          !< intermediate particle position (x)
151    REAL(wp) ::  pos_x_old      !< particle position (x) at previous particle timestep
152    REAL(wp) ::  pos_y          !< intermediate particle position (y)
153    REAL(wp) ::  pos_y_old      !< particle position (y) at previous particle timestep
154    REAL(wp) ::  pos_z          !< intermediate particle position (z)
155    REAL(wp) ::  pos_z_old      !< particle position (z) at previous particle timestep
156    REAL(wp) ::  prt_x          !< current particle position (x)
157    REAL(wp) ::  prt_y          !< current particle position (y)
158    REAL(wp) ::  prt_z          !< current particle position (z)
159    REAL(wp) ::  t_old          !< previous reflection time
160    REAL(wp) ::  tmp_t          !< dummy for sorting algorithmus
161    REAL(wp) ::  xwall          !< location of wall in x
162    REAL(wp) ::  ywall          !< location of wall in y
163    REAL(wp) ::  zwall1         !< location of wall in z (old grid box)
164    REAL(wp) ::  zwall2         !< location of wall in z (current grid box)
165    REAL(wp) ::  zwall3         !< location of wall in z (old y, current x)
166    REAL(wp) ::  zwall4         !< location of wall in z (current y, old x)
167
168    REAL(wp), DIMENSION(0:10) ::  t  !< reflection time
169
170
171    IF ( location == 'bottom/top' )  THEN
172
173!
174!--    Apply boundary conditions to those particles that have crossed the top or
175!--    bottom boundary and delete those particles, which are older than allowed
176       DO  n = 1, number_of_particles
177
178!
179!--       Stop if particles have moved further than the length of one
180!--       PE subdomain (newly released particles have age = age_m!)
181          IF ( particles(n)%age /= particles(n)%age_m )  THEN
182             IF ( ABS(particles(n)%speed_x) >                                  &
183                  ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m)  .OR. &
184                  ABS(particles(n)%speed_y) >                                  &
185                  ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) )  THEN
186
187                  WRITE( message_string, * )  'particle too fast.  n = ',  n
188                  CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
189             ENDIF
190          ENDIF
191
192          IF ( particles(n)%age > particle_maximum_age  .AND.  &
193               particles(n)%particle_mask )                              &
194          THEN
195             particles(n)%particle_mask  = .FALSE.
196             deleted_particles = deleted_particles + 1
197          ENDIF
198
199          IF ( particles(n)%z >= zu(nz)  .AND.  particles(n)%particle_mask )  THEN
200             IF ( ibc_par_t == 1 )  THEN
201!
202!--             Particle absorption
203                particles(n)%particle_mask  = .FALSE.
204                deleted_particles = deleted_particles + 1
205             ELSEIF ( ibc_par_t == 2 )  THEN
206!
207!--             Particle reflection
208                particles(n)%z       = 2.0_wp * zu(nz) - particles(n)%z
209                particles(n)%speed_z = -particles(n)%speed_z
210                IF ( use_sgs_for_particles  .AND. &
211                     particles(n)%rvar3 > 0.0_wp )  THEN
212                   particles(n)%rvar3 = -particles(n)%rvar3
213                ENDIF
214             ENDIF
215          ENDIF
216         
217          IF ( particles(n)%z < zw(0)  .AND.  particles(n)%particle_mask )  THEN
218             IF ( ibc_par_b == 1 )  THEN
219!
220!--             Particle absorption
221                particles(n)%particle_mask  = .FALSE.
222                deleted_particles = deleted_particles + 1
223             ELSEIF ( ibc_par_b == 2 )  THEN
224!
225!--             Particle reflection
226                particles(n)%z       = 2.0_wp * zw(0) - particles(n)%z
227                particles(n)%speed_z = -particles(n)%speed_z
228                IF ( use_sgs_for_particles  .AND. &
229                     particles(n)%rvar3 < 0.0_wp )  THEN
230                   particles(n)%rvar3 = -particles(n)%rvar3
231                ENDIF
232             ENDIF
233          ENDIF
234       ENDDO
235
236    ELSEIF ( location == 'walls' )  THEN
237
238
239       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
240
241       DO  n = 1, number_of_particles
242!
243!--       Recalculate particle timestep
244          dt_particle = particles(n)%age - particles(n)%age_m
245!
246!--       Obtain x/y indices for current particle position
247          i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
248          j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
249!
250!--       Save current particle positions
251          prt_x = particles(n)%x
252          prt_y = particles(n)%y
253          prt_z = particles(n)%z
254!
255!--       Recalculate old particle positions
256          pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
257          pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
258          pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
259!
260!--       Obtain x/y indices for old particle positions
261          i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
262          j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
263!
264!--       Determine horizontal as well as vertical walls at which particle can
265!--       be potentially reflected.
266!--       Start with walls aligned in yz layer.
267!--       Wall to the right
268          IF ( prt_x > pos_x_old )  THEN
269             xwall = ( i1 + 0.5_wp ) * dx
270!
271!--       Wall to the left
272          ELSE
273             xwall = ( i1 - 0.5_wp ) * dx
274          ENDIF
275!
276!--       Walls aligned in xz layer
277!--       Wall to the north
278          IF ( prt_y > pos_y_old )  THEN
279             ywall = ( j1 + 0.5_wp ) * dy
280!--       Wall to the south
281          ELSE
282             ywall = ( j1 - 0.5_wp ) * dy
283          ENDIF
284!
285!--       Walls aligned in xy layer at which particle can be possiblly reflected.
286!--       The construct of MERGE and BTEST is used to determine the topography-
287!--       top index (former nzb_s_inner).
288          zwall1 = zw( MAXLOC(                                                 &
289                          MERGE( 1, 0,                                         &
290                                 BTEST( wall_flags_0(nzb:nzb_max,j2,i2), 12 )  &
291                               ), DIM = 1                                      &
292                             ) - 1 )                                             
293          zwall2 = zw( MAXLOC(                                                 &
294                          MERGE( 1, 0,                                         &
295                                 BTEST( wall_flags_0(nzb:nzb_max,j1,i1), 12 )  &
296                               ), DIM = 1                                      &
297                             ) - 1 ) 
298          zwall3 = zw( MAXLOC(                                                 &
299                          MERGE( 1, 0,                                         &
300                                 BTEST( wall_flags_0(nzb:nzb_max,j1,i2), 12 )  &
301                               ), DIM = 1                                      &
302                             ) - 1 ) 
303          zwall4 = zw( MAXLOC(                                                 &
304                          MERGE( 1, 0,                                         &
305                                 BTEST( wall_flags_0(nzb:nzb_max,j2,i1), 12 )  &
306                               ), DIM = 1                                      &
307                             ) - 1 ) 
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 = MAXLOC(                                               &
493                            MERGE( 1, 0,                                       &
494                                 BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 )  &
495                                 ), DIM = 1                                    &
496                               ) - 1 
497                IF ( ABS( pos_x - xwall ) < eps      .AND.                     &
498                     pos_z <= zw(k_wall)             .AND.                     &
499                     reach_x(t_index)                .AND.                     &
500                     .NOT. reflect_x )  THEN
501!
502!--                Reflection in x-direction.
503!--                Ensure correct reflection by MIN/MAX functions, depending on
504!--                direction of particle transport.
505!--                Due to rounding errors pos_x do not exactly matches the wall
506!--                location, leading to erroneous reflection.             
507                   pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ),        &
508                                  MAX( 2.0_wp * xwall - pos_x, xwall ),        &
509                                  particles(n)%x > xwall )
510!
511!--                Change sign of particle speed                     
512                   particles(n)%speed_x = - particles(n)%speed_x
513!
514!--                Change also sign of subgrid-scale particle speed
515                   particles(n)%rvar1 = - particles(n)%rvar1
516!
517!--                Set flag that reflection along x is already done
518                   reflect_x          = .TRUE.
519!
520!--                As particle do not crosses any further yz-wall during
521!--                this timestep, set further x-indices to the current one.
522                   x_ind(t_index:t_index_number) = i1
523!
524!--             If particle already reached the wall but was not reflected,
525!--             set further x-indices to the new one.
526                ELSEIF ( x_wall_reached .AND. .NOT. reflect_x )  THEN
527                    x_ind(t_index:t_index_number) = i2
528                ENDIF
529!
530!--             Check if a particle needs to be reflected at any xz-wall. If
531!--             necessary, carry out reflection. At first, determine index of
532!--             topography top at (j3,i3) 
533                k_wall = MAXLOC(                                               &
534                            MERGE( 1, 0,                                       &
535                                 BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 )  &
536                                 ), DIM = 1                                    &
537                               ) - 1 
538                IF ( ABS( pos_y - ywall ) < eps      .AND.                     &
539                     pos_z <= zw(k_wall)             .AND.                     &
540                     reach_y(t_index)                .AND.                     &
541                     .NOT. reflect_y ) THEN
542
543                   pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ),        &
544                                  MAX( 2.0_wp * ywall - pos_y, ywall ),        &
545                                  particles(n)%y > ywall ) 
546
547                   particles(n)%speed_y = - particles(n)%speed_y     
548                   particles(n)%rvar2   = - particles(n)%rvar2       
549     
550                   reflect_y            = .TRUE.
551                   y_ind(t_index:t_index_number) = j1
552
553                ELSEIF ( y_wall_reached .AND. .NOT. reflect_y )  THEN
554                   y_ind(t_index:t_index_number) = j2
555                ENDIF
556!
557!--             Check if a particle needs to be reflected at any xy-wall. If
558!--             necessary, carry out reflection.
559                IF ( downwards .AND. reach_z(t_index) .AND.                    &
560                     .NOT. reflect_z )  THEN
561!
562!--                Determine index of topography top at (j3,i3) and chick if
563!--                particle is below. 
564                   k_wall = MAXLOC(                                            &
565                          MERGE( 1, 0,                                         &
566                                   BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 )&
567                               ), DIM = 1                                      &
568                                  ) - 1 
569                   IF ( pos_z - zw(k_wall) < eps ) THEN
570 
571                      pos_z = MAX( 2.0_wp * zw(k_wall) - pos_z,    &
572                                   zw(k_wall) )
573
574                      particles(n)%speed_z = - particles(n)%speed_z
575                      particles(n)%rvar3   = - particles(n)%rvar3
576
577                      reflect_z            = .TRUE.
578
579                   ENDIF
580
581                ENDIF
582!
583!--             Swap time
584                t_old = t(t_index)
585
586             ENDDO
587!
588!--          If a particle was reflected, calculate final position from last
589!--          intermediate position.
590             IF ( reflect_x .OR. reflect_y .OR. reflect_z )  THEN
591
592                particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle      &
593                                                         * particles(n)%speed_x
594                particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle      &
595                                                         * particles(n)%speed_y
596                particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle      &
597                                                         * particles(n)%speed_z
598
599             ENDIF
600
601          ENDIF
602
603       ENDDO
604
605       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
606
607    ENDIF
608
609 END SUBROUTINE lpm_boundary_conds
Note: See TracBrowser for help on using the repository browser.