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

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

get topograpyh top index via function call

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