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

Last change on this file since 858 was 850, checked in by raasch, 13 years ago

last commit documented

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