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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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