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

Last change on this file since 2258 was 2233, checked in by suehring, 7 years ago

last commit documented

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