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

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

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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