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

Last change on this file since 1319 was 1319, checked in by raasch, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 24.7 KB
Line 
1 SUBROUTINE lpm_boundary_conds( range )
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_boundary_conds.f90 1319 2014-03-17 15:08:44Z raasch $
27!
28! 1318 2014-03-17 13:35:16Z raasch
29! module interfaces removed
30!
31! 1036 2012-10-22 13:43:42Z raasch
32! code put under GPL (PALM 3.9)
33!
34! 849 2012-03-15 10:35:09Z raasch
35! routine renamed lpm_boundary_conds, bottom and top boundary conditions
36! included (former part of advec_particles)
37!
38! 824 2012-02-17 09:09:57Z raasch
39! particle attributes speed_x|y|z_sgs renamed rvar1|2|3
40!
41! 150 2008-02-29 08:19:58Z raasch
42! Vertical index calculations adjusted for ocean runs.
43!
44! Initial version (2007/03/09)
45!
46! Description:
47! ------------
48! Boundary conditions for the Lagrangian particles.
49! The routine consists of two different parts. One handles the bottom (flat)
50! and top boundary. In this part, also particles which exceeded their lifetime
51! are deleted.
52! The other part handles the reflection of particles from vertical walls.
53! This part was developed by Jin Zhang during 2006-2007.
54!
55! To do: Code structure for finding the t_index values and for checking the
56! -----  reflection conditions is basically the same for all four cases, so it
57!        should be possible to further simplify/shorten it.
58!
59! THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
60! (see offset_ocean_*)
61!------------------------------------------------------------------------------!
62
63    USE arrays_3d
64    USE control_parameters
65    USE cpulog
66    USE grid_variables
67    USE indices
68    USE particle_attributes
69    USE pegrid
70
71    IMPLICIT NONE
72
73    CHARACTER (LEN=*) ::  range
74
75    INTEGER ::  i, inc, ir, i1, i2, i3, i5, j, jr, j1, j2, j3, j5, k, k1, k2, &
76                k3, k5, n, nn, t_index, t_index_number
77
78    LOGICAL ::  reflect_x, reflect_y, reflect_z
79
80    REAL ::  dt_particle, pos_x, pos_x_old, pos_y, pos_y_old, pos_z, &
81             pos_z_old, prt_x, prt_y, prt_z, tmp_t, xline, yline, zline
82
83    REAL ::  t(1:200)
84
85
86
87    IF ( range == 'bottom/top' )  THEN
88
89!
90!--    Apply boundary conditions to those particles that have crossed the top or
91!--    bottom boundary and delete those particles, which are older than allowed
92       DO  n = 1, number_of_particles
93
94          nn = particles(n)%tail_id
95
96!
97!--       Stop if particles have moved further than the length of one
98!--       PE subdomain (newly released particles have age = age_m!)
99          IF ( particles(n)%age /= particles(n)%age_m )  THEN
100             IF ( ABS(particles(n)%speed_x) >                                  &
101                  ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m)  .OR. &
102                  ABS(particles(n)%speed_y) >                                  &
103                  ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) )  THEN
104
105                  WRITE( message_string, * )  'particle too fast.  n = ',  n
106                  CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
107             ENDIF
108          ENDIF
109
110          IF ( particles(n)%age > particle_maximum_age  .AND.  &
111               particle_mask(n) )                              &
112          THEN
113             particle_mask(n)  = .FALSE.
114             deleted_particles = deleted_particles + 1
115             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
116                tail_mask(nn) = .FALSE.
117                deleted_tails = deleted_tails + 1
118             ENDIF
119          ENDIF
120
121          IF ( particles(n)%z >= zu(nz)  .AND.  particle_mask(n) )  THEN
122             IF ( ibc_par_t == 1 )  THEN
123!
124!--             Particle absorption
125                particle_mask(n)  = .FALSE.
126                deleted_particles = deleted_particles + 1
127                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
128                   tail_mask(nn) = .FALSE.
129                   deleted_tails = deleted_tails + 1
130                ENDIF
131             ELSEIF ( ibc_par_t == 2 )  THEN
132!
133!--             Particle reflection
134                particles(n)%z       = 2.0 * zu(nz) - particles(n)%z
135                particles(n)%speed_z = -particles(n)%speed_z
136                IF ( use_sgs_for_particles  .AND. &
137                     particles(n)%rvar3 > 0.0 )  THEN
138                   particles(n)%rvar3 = -particles(n)%rvar3
139                ENDIF
140                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
141                   particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
142                                               particle_tail_coordinates(1,3,nn)
143                ENDIF
144             ENDIF
145          ENDIF
146          IF ( particles(n)%z < zw(0)  .AND.  particle_mask(n) )  THEN
147             IF ( ibc_par_b == 1 )  THEN
148!
149!--             Particle absorption
150                particle_mask(n)  = .FALSE.
151                deleted_particles = deleted_particles + 1
152                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
153                   tail_mask(nn) = .FALSE.
154                   deleted_tails = deleted_tails + 1
155                ENDIF
156             ELSEIF ( ibc_par_b == 2 )  THEN
157!
158!--             Particle reflection
159                particles(n)%z       = 2.0 * zw(0) - particles(n)%z
160                particles(n)%speed_z = -particles(n)%speed_z
161                IF ( use_sgs_for_particles  .AND. &
162                     particles(n)%rvar3 < 0.0 )  THEN
163                   particles(n)%rvar3 = -particles(n)%rvar3
164                ENDIF
165                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
166                   particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
167                                               particle_tail_coordinates(1,3,nn)
168                ENDIF
169                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
170                   particle_tail_coordinates(1,3,nn) = 2.0 * zw(0) - &
171                                               particle_tail_coordinates(1,3,nn)
172                ENDIF
173             ENDIF
174          ENDIF
175       ENDDO
176
177    ELSEIF ( range == 'walls' )  THEN
178
179       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
180
181       reflect_x = .FALSE.
182       reflect_y = .FALSE.
183       reflect_z = .FALSE.
184
185       DO  n = 1, number_of_particles
186
187          dt_particle = particles(n)%age - particles(n)%age_m
188
189          i2 = ( particles(n)%x + 0.5 * dx ) * ddx
190          j2 = ( particles(n)%y + 0.5 * dy ) * ddy
191          k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
192
193          prt_x = particles(n)%x
194          prt_y = particles(n)%y
195          prt_z = particles(n)%z
196
197!
198!--       If the particle position is below the surface, it has to be reflected
199          IF ( k2 <= nzb_s_inner(j2,i2)  .AND.  nzb_s_inner(j2,i2) /=0 )  THEN
200
201             pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
202             pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
203             pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
204             i1 = ( pos_x_old + 0.5 * dx ) * ddx
205             j1 = ( pos_y_old + 0.5 * dy ) * ddy
206             k1 = pos_z_old / dz + offset_ocean_nzt_m1
207
208!
209!--          Case 1
210             IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old )&
211             THEN
212                t_index = 1
213
214                DO  i = i1, i2
215                   xline      = i * dx + 0.5 * dx
216                   t(t_index) = ( xline - pos_x_old ) / &
217                                ( particles(n)%x - pos_x_old )
218                   t_index    = t_index + 1
219                ENDDO
220
221                DO  j = j1, j2
222                   yline      = j * dy + 0.5 * dy
223                   t(t_index) = ( yline - pos_y_old ) / &
224                                ( particles(n)%y - pos_y_old )
225                   t_index    = t_index + 1
226                ENDDO
227
228                IF ( particles(n)%z < pos_z_old )  THEN
229                   DO  k = k1, k2 , -1
230                      zline      = k * dz
231                      t(t_index) = ( pos_z_old - zline ) / &
232                                   ( pos_z_old - particles(n)%z )
233                      t_index    = t_index + 1
234                   ENDDO
235                ENDIF
236
237                t_index_number = t_index - 1
238
239!
240!--             Sorting t
241                inc = 1
242                jr  = 1
243                DO WHILE ( inc <= t_index_number )
244                   inc = 3 * inc + 1
245                ENDDO
246
247                DO WHILE ( inc > 1 )
248                   inc = inc / 3
249                   DO  ir = inc+1, t_index_number
250                      tmp_t = t(ir)
251                      jr    = ir
252                      DO WHILE ( t(jr-inc) > tmp_t )
253                         t(jr) = t(jr-inc)
254                         jr    = jr - inc
255                         IF ( jr <= inc )  EXIT
256                      ENDDO
257                      t(jr) = tmp_t
258                   ENDDO
259                ENDDO
260
261         case1: DO  t_index = 1, t_index_number
262
263                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
264                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
265                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
266
267                   i3 = ( pos_x + 0.5 * dx ) * ddx   
268                   j3 = ( pos_y + 0.5 * dy ) * ddy
269                   k3 = pos_z / dz + offset_ocean_nzt_m1
270
271                   i5 = pos_x * ddx
272                   j5 = pos_y * ddy
273                   k5 = pos_z / dz + offset_ocean_nzt_m1
274
275                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
276                        nzb_s_inner(j5,i3) /= 0 )  THEN
277
278                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
279                           k3 == nzb_s_inner(j5,i3) )  THEN
280                         reflect_z = .TRUE.
281                         EXIT case1
282                      ENDIF
283
284                   ENDIF
285
286                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
287                        nzb_s_inner(j3,i5) /= 0 )  THEN
288
289                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
290                           k3 == nzb_s_inner(j3,i5) )  THEN
291                         reflect_z = .TRUE.
292                         EXIT case1
293                      ENDIF
294
295                   ENDIF
296
297                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
298                        nzb_s_inner(j3,i3) /= 0 )  THEN
299
300                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
301                           k3 == nzb_s_inner(j3,i3) )  THEN
302                         reflect_z = .TRUE.
303                         EXIT case1
304                      ENDIF
305
306                      IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
307                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
308                         reflect_y = .TRUE.
309                         EXIT case1
310                      ENDIF
311
312                      IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
313                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
314                         reflect_x = .TRUE.
315                         EXIT case1
316                      ENDIF
317
318                   ENDIF
319
320                ENDDO case1
321
322!
323!--          Case 2
324             ELSEIF ( particles(n)%x > pos_x_old  .AND. &
325                      particles(n)%y < pos_y_old )  THEN
326
327                t_index = 1
328
329                DO  i = i1, i2
330                   xline      = i * dx + 0.5 * dx
331                   t(t_index) = ( xline - pos_x_old ) / &
332                                ( particles(n)%x - pos_x_old )
333                   t_index    = t_index + 1
334                ENDDO
335
336                DO  j = j1, j2, -1
337                   yline      = j * dy - 0.5 * dy
338                   t(t_index) = ( pos_y_old - yline ) / &
339                                ( pos_y_old - particles(n)%y )
340                   t_index    = t_index + 1
341                ENDDO
342
343                IF ( particles(n)%z < pos_z_old )  THEN
344                   DO  k = k1, k2 , -1
345                      zline      = k * dz
346                      t(t_index) = ( pos_z_old - zline ) / &
347                                   ( pos_z_old - particles(n)%z )
348                      t_index    = t_index + 1
349                   ENDDO
350                ENDIF
351                t_index_number = t_index-1
352
353!
354!--             Sorting t
355                inc = 1
356                jr  = 1
357                DO WHILE ( inc <= t_index_number )
358                   inc = 3 * inc + 1
359                ENDDO
360
361                DO WHILE ( inc > 1 )
362                   inc = inc / 3
363                   DO  ir = inc+1, t_index_number
364                      tmp_t = t(ir)
365                      jr    = ir
366                      DO WHILE ( t(jr-inc) > tmp_t )
367                         t(jr) = t(jr-inc)
368                         jr    = jr - inc
369                         IF ( jr <= inc )  EXIT
370                      ENDDO
371                      t(jr) = tmp_t
372                   ENDDO
373                ENDDO
374
375         case2: DO  t_index = 1, t_index_number
376
377                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
378                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
379                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
380
381                   i3 = ( pos_x + 0.5 * dx ) * ddx
382                   j3 = ( pos_y + 0.5 * dy ) * ddy
383                   k3 = pos_z / dz + offset_ocean_nzt_m1
384
385                   i5 = pos_x * ddx
386                   j5 = pos_y * ddy
387                   k5 = pos_z / dz + offset_ocean_nzt_m1
388
389                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
390                        nzb_s_inner(j3,i5) /= 0 )  THEN
391
392                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
393                           k3 == nzb_s_inner(j3,i5) )  THEN
394                         reflect_z = .TRUE.
395                         EXIT case2
396                      ENDIF
397
398                   ENDIF
399
400                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
401                        nzb_s_inner(j3,i3) /= 0 )  THEN
402
403                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
404                           k3 == nzb_s_inner(j3,i3) )  THEN
405                         reflect_z = .TRUE.
406                         EXIT case2
407                      ENDIF
408
409                      IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
410                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
411                         reflect_x = .TRUE.
412                         EXIT case2
413                      ENDIF
414
415                   ENDIF
416
417                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
418                        nzb_s_inner(j5,i3) /= 0 )  THEN
419
420                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
421                           k3 == nzb_s_inner(j5,i3) )  THEN
422                         reflect_z = .TRUE.
423                         EXIT case2
424                      ENDIF
425
426                      IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
427                           pos_z < nzb_s_inner(j5,i3) * dz )  THEN
428                         reflect_y = .TRUE.
429                         EXIT case2
430                      ENDIF
431
432                   ENDIF
433
434                ENDDO case2
435
436!
437!--          Case 3
438             ELSEIF ( particles(n)%x < pos_x_old  .AND. &
439                      particles(n)%y > pos_y_old )  THEN
440
441                t_index = 1
442
443                DO  i = i1, i2, -1
444                   xline      = i * dx - 0.5 * dx
445                   t(t_index) = ( pos_x_old - xline ) / &
446                                ( pos_x_old - particles(n)%x )
447                   t_index    = t_index + 1
448                ENDDO
449
450                DO  j = j1, j2
451                   yline      = j * dy + 0.5 * dy
452                   t(t_index) = ( yline - pos_y_old ) / &
453                                ( particles(n)%y - pos_y_old )
454                   t_index    = t_index + 1
455                ENDDO
456
457                IF ( particles(n)%z < pos_z_old )  THEN
458                   DO  k = k1, k2 , -1
459                      zline      = k * dz
460                      t(t_index) = ( pos_z_old - zline ) / &
461                                   ( pos_z_old - particles(n)%z )
462                      t_index    = t_index + 1
463                   ENDDO
464                ENDIF
465                t_index_number = t_index - 1
466
467!
468!--             Sorting t
469                inc = 1
470                jr  = 1
471
472                DO WHILE ( inc <= t_index_number )
473                   inc = 3 * inc + 1
474                ENDDO
475
476                DO WHILE ( inc > 1 )
477                   inc = inc / 3
478                   DO  ir = inc+1, t_index_number
479                      tmp_t = t(ir)
480                      jr    = ir
481                      DO WHILE ( t(jr-inc) > tmp_t )
482                         t(jr) = t(jr-inc)
483                         jr    = jr - inc
484                         IF ( jr <= inc )  EXIT
485                      ENDDO
486                      t(jr) = tmp_t
487                   ENDDO
488                ENDDO
489
490         case3: DO  t_index = 1, t_index_number
491
492                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
493                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
494                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
495
496                   i3 = ( pos_x + 0.5 * dx ) * ddx
497                   j3 = ( pos_y + 0.5 * dy ) * ddy
498                   k3 = pos_z / dz + offset_ocean_nzt_m1
499
500                   i5 = pos_x * ddx
501                   j5 = pos_y * ddy
502                   k5 = pos_z / dz + offset_ocean_nzt_m1
503
504                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
505                        nzb_s_inner(j5,i3) /= 0 )  THEN
506
507                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
508                           k3 == nzb_s_inner(j5,i3) )  THEN
509                         reflect_z = .TRUE.
510                         EXIT case3
511                      ENDIF
512
513                   ENDIF
514
515                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
516                        nzb_s_inner(j3,i3) /= 0 )  THEN
517
518                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
519                           k3 == nzb_s_inner(j3,i3) )  THEN
520                         reflect_z = .TRUE.
521                         EXIT case3
522                      ENDIF
523
524                      IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
525                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
526                         reflect_y = .TRUE.
527                         EXIT case3
528                      ENDIF
529
530                   ENDIF
531
532                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
533                        nzb_s_inner(j3,i5) /= 0 )  THEN
534
535                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
536                           k3 == nzb_s_inner(j3,i5) )  THEN
537                         reflect_z = .TRUE.
538                         EXIT case3
539                      ENDIF
540
541                      IF ( pos_x == ( i5 * dx + 0.5 * dx )  .AND. &
542                           pos_z < nzb_s_inner(j3,i5) * dz )  THEN
543                         reflect_x = .TRUE.
544                         EXIT case3
545                      ENDIF
546
547                   ENDIF
548
549                ENDDO case3
550
551!
552!--          Case 4
553             ELSEIF ( particles(n)%x < pos_x_old  .AND. &
554                      particles(n)%y < pos_y_old )  THEN
555
556                t_index = 1
557
558                DO  i = i1, i2, -1
559                   xline      = i * dx - 0.5 * dx
560                   t(t_index) = ( pos_x_old - xline ) / &
561                                ( pos_x_old - particles(n)%x )
562                   t_index    = t_index + 1
563                ENDDO
564
565                DO  j = j1, j2, -1
566                   yline      = j * dy - 0.5 * dy
567                   t(t_index) = ( pos_y_old - yline ) / &
568                                ( pos_y_old - particles(n)%y )
569                   t_index    = t_index + 1
570                ENDDO
571
572                IF ( particles(n)%z < pos_z_old )  THEN
573                   DO  k = k1, k2 , -1
574                      zline      = k * dz
575                      t(t_index) = ( pos_z_old - zline ) / &
576                                   ( pos_z_old-particles(n)%z )
577                      t_index    = t_index + 1
578                   ENDDO
579                ENDIF
580                t_index_number = t_index-1
581
582!
583!--             Sorting t
584                inc = 1
585                jr  = 1
586
587                DO WHILE ( inc <= t_index_number )
588                   inc = 3 * inc + 1
589                ENDDO
590
591                DO WHILE ( inc > 1 )
592                   inc = inc / 3
593                   DO  ir = inc+1, t_index_number
594                      tmp_t = t(ir)
595                      jr = ir
596                      DO WHILE ( t(jr-inc) > tmp_t )
597                         t(jr) = t(jr-inc)
598                         jr    = jr - inc
599                         IF ( jr <= inc )  EXIT
600                      ENDDO
601                      t(jr) = tmp_t
602                   ENDDO
603                ENDDO
604
605         case4: DO  t_index = 1, t_index_number
606
607                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
608                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
609                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
610
611                   i3 = ( pos_x + 0.5 * dx ) * ddx   
612                   j3 = ( pos_y + 0.5 * dy ) * ddy
613                   k3 = pos_z / dz + offset_ocean_nzt_m1
614
615                   i5 = pos_x * ddx
616                   j5 = pos_y * ddy
617                   k5 = pos_z / dz + offset_ocean_nzt_m1
618
619                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
620                        nzb_s_inner(j3,i3) /= 0 )  THEN
621
622                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
623                           k3 == nzb_s_inner(j3,i3) )  THEN
624                         reflect_z = .TRUE.
625                         EXIT case4
626                      ENDIF
627
628                   ENDIF
629
630                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
631                        nzb_s_inner(j3,i5) /= 0 )  THEN
632
633                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
634                           k3 == nzb_s_inner(j3,i5) )  THEN
635                         reflect_z = .TRUE.
636                         EXIT case4
637                      ENDIF
638
639                      IF ( pos_x == ( i5 * dx + 0.5 * dx )  .AND. &
640                           nzb_s_inner(j3,i5) /=0  .AND.          &
641                           pos_z < nzb_s_inner(j3,i5) * dz )  THEN
642                         reflect_x = .TRUE.
643                         EXIT case4
644                      ENDIF
645
646                   ENDIF
647
648                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
649                        nzb_s_inner(j5,i3) /= 0 )  THEN
650
651                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
652                           k5 == nzb_s_inner(j5,i3) )  THEN
653                         reflect_z = .TRUE.
654                         EXIT case4
655                      ENDIF
656
657                      IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
658                           nzb_s_inner(j5,i3) /= 0  .AND.         &
659                           pos_z < nzb_s_inner(j5,i3) * dz )  THEN
660                         reflect_y = .TRUE.
661                         EXIT case4
662                      ENDIF
663
664                   ENDIF
665 
666                ENDDO case4
667
668             ENDIF
669
670          ENDIF      ! Check, if particle position is below the surface
671
672!
673!--       Do the respective reflection, in case that one of the above conditions
674!--       is found to be true
675          IF ( reflect_z )  THEN
676
677             particles(n)%z       = 2.0 * pos_z - prt_z
678             particles(n)%speed_z = - particles(n)%speed_z
679
680             IF ( use_sgs_for_particles )  THEN
681                particles(n)%rvar3 = - particles(n)%rvar3
682             ENDIF
683             reflect_z = .FALSE.
684
685          ELSEIF ( reflect_y )  THEN
686
687             particles(n)%y       = 2.0 * pos_y - prt_y
688             particles(n)%speed_y = - particles(n)%speed_y
689
690             IF ( use_sgs_for_particles )  THEN
691                particles(n)%rvar2 = - particles(n)%rvar2
692             ENDIF
693             reflect_y = .FALSE.
694
695          ELSEIF ( reflect_x )  THEN
696
697             particles(n)%x       = 2.0 * pos_x - prt_x
698             particles(n)%speed_x = - particles(n)%speed_x
699
700             IF ( use_sgs_for_particles )  THEN
701                particles(n)%rvar1 = - particles(n)%rvar1
702             ENDIF
703             reflect_x = .FALSE.
704
705          ENDIF       
706   
707       ENDDO
708
709       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
710
711    ENDIF
712
713 END SUBROUTINE lpm_boundary_conds
Note: See TracBrowser for help on using the repository browser.