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

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

last commit documented

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