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

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

last commit documented

  • Property svn:keywords set to Id
File size: 23.8 KB
RevLine 
[849]1 SUBROUTINE lpm_boundary_conds( range )
[58]2
3!------------------------------------------------------------------------------!
[484]4! Current revisions:
[58]5! -----------------
[850]6!
[58]7!
8! Former revisions:
9! -----------------
10! $Id: lpm_boundary_conds.f90 850 2012-03-15 12:09:25Z raasch $
[198]11!
[850]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!
[826]16! 824 2012-02-17 09:09:57Z raasch
17! particle attributes speed_x|y|z_sgs renamed rvar1|2|3
18!
[198]19! 150 2008-02-29 08:19:58Z raasch
20! Vertical index calculations adjusted for ocean runs.
21!
[58]22! Initial version (2007/03/09)
23!
24! Description:
25! ------------
[849]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!
[61]33! To do: Code structure for finding the t_index values and for checking the
[849]34! -----  reflection conditions is basically the same for all four cases, so it
[61]35!        should be possible to further simplify/shorten it.
[849]36!
37! THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
38! (see offset_ocean_*)
[58]39!------------------------------------------------------------------------------!
40
[849]41    USE arrays_3d
[58]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
[849]52    CHARACTER (LEN=*) ::  range
53
[60]54    INTEGER ::  i, inc, ir, i1, i2, i3, i5, j, jr, j1, j2, j3, j5, k, k1, k2, &
[849]55                k3, k5, n, nn, t_index, t_index_number
[58]56
[61]57    LOGICAL ::  reflect_x, reflect_y, reflect_z
58
[60]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
[58]61
[60]62    REAL ::  t(1:200)
63
[58]64
65
[849]66    IF ( range == 'bottom/top' )  THEN
[58]67
[849]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
[61]72
[849]73          nn = particles(n)%tail_id
[58]74
[849]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
[60]83
[849]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
[58]88
[849]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
[58]99
[849]100          IF ( particles(n)%z >= zu(nz)  .AND.  particle_mask(n) )  THEN
101             IF ( ibc_par_t == 1 )  THEN
[61]102!
[849]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
[58]155
[849]156    ELSEIF ( range == 'walls' )  THEN
[58]157
[849]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
[58]176!
[849]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
[58]179
[849]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
[58]186
[849]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
[58]192
[849]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 )
[58]197                   t_index    = t_index + 1
198                ENDDO
199
[849]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
[58]206
[849]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
[58]218!
[849]219!--             Sorting t
220                inc = 1
221                jr  = 1
222                DO WHILE ( inc <= t_index_number )
223                   inc = 3 * inc + 1
224                ENDDO
[58]225
[849]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
[58]237                   ENDDO
238                ENDDO
239
[849]240         case1: DO  t_index = 1, t_index_number
[58]241
[849]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 )
[58]245
[849]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
[58]249
[849]250                   i5 = pos_x * ddx
251                   j5 = pos_y * ddy
252                   k5 = pos_z / dz + offset_ocean_nzt_m1
[58]253
[849]254                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
255                        nzb_s_inner(j5,i3) /= 0 )  THEN
[61]256
[849]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
[61]263                   ENDIF
[58]264
[849]265                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
266                        nzb_s_inner(j3,i5) /= 0 )  THEN
[58]267
[849]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
[61]273
274                   ENDIF
[58]275
[849]276                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
277                        nzb_s_inner(j3,i3) /= 0 )  THEN
[58]278
[849]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
[61]284
[849]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
[58]290
[849]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
[58]296
[61]297                   ENDIF
[58]298
[849]299                ENDDO case1
[58]300
301!
[849]302!--          Case 2
303             ELSEIF ( particles(n)%x > pos_x_old  .AND. &
304                      particles(n)%y < pos_y_old )  THEN
[61]305
[849]306                t_index = 1
[58]307
[849]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
[58]314
[849]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 )
[58]319                   t_index    = t_index + 1
320                ENDDO
321
[849]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
[58]332!
[849]333!--             Sorting t
334                inc = 1
335                jr  = 1
336                DO WHILE ( inc <= t_index_number )
337                   inc = 3 * inc + 1
338                ENDDO
[58]339
[849]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
[58]351                   ENDDO
352                ENDDO
353
[849]354         case2: DO  t_index = 1, t_index_number
[58]355
[849]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 )
[58]359
[849]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
[58]363
[849]364                   i5 = pos_x * ddx
365                   j5 = pos_y * ddy
366                   k5 = pos_z / dz + offset_ocean_nzt_m1
[58]367
[849]368                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
369                        nzb_s_inner(j3,i5) /= 0 )  THEN
[61]370
[849]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
[61]377                   ENDIF
[58]378
[849]379                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
380                        nzb_s_inner(j3,i3) /= 0 )  THEN
[58]381
[849]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
[61]387
[849]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
[58]393
[61]394                   ENDIF
[58]395
[849]396                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
397                        nzb_s_inner(j5,i3) /= 0 )  THEN
[58]398
[849]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
[61]404
[849]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
[58]410
[61]411                   ENDIF
[58]412
[849]413                ENDDO case2
[58]414
415!
[849]416!--          Case 3
417             ELSEIF ( particles(n)%x < pos_x_old  .AND. &
418                      particles(n)%y > pos_y_old )  THEN
[61]419
[849]420                t_index = 1
[58]421
[849]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
[58]428
[849]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 )
[58]433                   t_index    = t_index + 1
434                ENDDO
435
[849]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
[58]446!
[849]447!--             Sorting t
448                inc = 1
449                jr  = 1
[58]450
[849]451                DO WHILE ( inc <= t_index_number )
452                   inc = 3 * inc + 1
453                ENDDO
[58]454
[849]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
[58]466                   ENDDO
467                ENDDO
468
[849]469         case3: DO  t_index = 1, t_index_number
[58]470
[849]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 )
[58]474
[849]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
[58]478
[849]479                   i5 = pos_x * ddx
480                   j5 = pos_y * ddy
481                   k5 = pos_z / dz + offset_ocean_nzt_m1
[58]482
[849]483                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
484                        nzb_s_inner(j5,i3) /= 0 )  THEN
[61]485
[849]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
[61]492                   ENDIF
[58]493
[849]494                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
495                        nzb_s_inner(j3,i3) /= 0 )  THEN
[58]496
[849]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
[61]502
[849]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
[58]508
[61]509                   ENDIF
[58]510
[849]511                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
512                        nzb_s_inner(j3,i5) /= 0 )  THEN
[58]513
[849]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
[61]519
[849]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
[58]525
526                   ENDIF
527
[849]528                ENDDO case3
[61]529
[58]530!
[849]531!--          Case 4
532             ELSEIF ( particles(n)%x < pos_x_old  .AND. &
533                      particles(n)%y < pos_y_old )  THEN
[61]534
[849]535                t_index = 1
[58]536
[849]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
[58]543
[849]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 )
[58]548                   t_index    = t_index + 1
549                ENDDO
550
[849]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
[58]561!
[849]562!--             Sorting t
563                inc = 1
564                jr  = 1
[58]565
[849]566                DO WHILE ( inc <= t_index_number )
567                   inc = 3 * inc + 1
568                ENDDO
[58]569
[849]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
[58]581                   ENDDO
582                ENDDO
583
[849]584         case4: DO  t_index = 1, t_index_number
[58]585
[849]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 )
[58]589
[849]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
[58]593
[849]594                   i5 = pos_x * ddx
595                   j5 = pos_y * ddy
596                   k5 = pos_z / dz + offset_ocean_nzt_m1
[58]597
[849]598                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
599                        nzb_s_inner(j3,i3) /= 0 )  THEN
[61]600
[849]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
[61]607                   ENDIF
[58]608
[849]609                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
610                        nzb_s_inner(j3,i5) /= 0 )  THEN
[58]611
[849]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
[61]617
[849]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
[58]624
[61]625                   ENDIF
[58]626
[849]627                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
628                        nzb_s_inner(j5,i3) /= 0 )  THEN
[61]629
[849]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
[61]635
[849]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
[58]642
[61]643                   ENDIF
644 
[849]645                ENDDO case4
[61]646
[849]647             ENDIF
[58]648
[849]649          ENDIF      ! Check, if particle position is below the surface
[61]650
651!
[849]652!--       Do the respective reflection, in case that one of the above conditions
653!--       is found to be true
654          IF ( reflect_z )  THEN
[61]655
[849]656             particles(n)%z       = 2.0 * pos_z - prt_z
657             particles(n)%speed_z = - particles(n)%speed_z
[61]658
[849]659             IF ( use_sgs_for_particles )  THEN
660                particles(n)%rvar3 = - particles(n)%rvar3
661             ENDIF
662             reflect_z = .FALSE.
[61]663
[849]664          ELSEIF ( reflect_y )  THEN
[61]665
[849]666             particles(n)%y       = 2.0 * pos_y - prt_y
667             particles(n)%speed_y = - particles(n)%speed_y
[61]668
[849]669             IF ( use_sgs_for_particles )  THEN
670                particles(n)%rvar2 = - particles(n)%rvar2
671             ENDIF
672             reflect_y = .FALSE.
[61]673
[849]674          ELSEIF ( reflect_x )  THEN
[61]675
[849]676             particles(n)%x       = 2.0 * pos_x - prt_x
677             particles(n)%speed_x = - particles(n)%speed_x
[61]678
[849]679             IF ( use_sgs_for_particles )  THEN
680                particles(n)%rvar1 = - particles(n)%rvar1
681             ENDIF
682             reflect_x = .FALSE.
[61]683
[849]684          ENDIF       
[58]685   
[849]686       ENDDO
[58]687
[849]688       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
[58]689
[849]690    ENDIF
[61]691
[849]692 END SUBROUTINE lpm_boundary_conds
Note: See TracBrowser for help on using the repository browser.