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

Last change on this file since 1822 was 1822, checked in by hoffmann, 5 years ago

changes in LPM and bulk cloud microphysics

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