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

Last change on this file since 1585 was 1360, checked in by hoffmann, 10 years ago

last commit documented

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