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

Last change on this file since 1734 was 1683, checked in by knoop, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 27.2 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-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: lpm_boundary_conds.f90 1683 2015-10-07 23:57:51Z raasch $
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, simulated_time
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, deleted_tails, ibc_par_b, ibc_par_t,         &
92               number_of_particles, particles,                                 &
93               particle_tail_coordinates, particle_type, offset_ocean_nzt_m1,  &
94               tail_mask, use_particle_tails, use_sgs_for_particles
95
96    USE pegrid
97
98    IMPLICIT NONE
99
100    CHARACTER (LEN=*) ::  range     !<
101   
102    INTEGER(iwp) ::  i              !<
103    INTEGER(iwp) ::  inc            !<
104    INTEGER(iwp) ::  ir             !<
105    INTEGER(iwp) ::  i1             !<
106    INTEGER(iwp) ::  i2             !<
107    INTEGER(iwp) ::  i3             !<
108    INTEGER(iwp) ::  i5             !<
109    INTEGER(iwp) ::  j              !<
110    INTEGER(iwp) ::  jr             !<
111    INTEGER(iwp) ::  j1             !<
112    INTEGER(iwp) ::  j2             !<
113    INTEGER(iwp) ::  j3             !<
114    INTEGER(iwp) ::  j5             !<
115    INTEGER(iwp) ::  k              !<
116    INTEGER(iwp) ::  k1             !<
117    INTEGER(iwp) ::  k2             !<
118    INTEGER(iwp) ::  k3             !<
119    INTEGER(iwp) ::  k5             !<
120    INTEGER(iwp) ::  n              !<
121    INTEGER(iwp) ::  nn             !<
122    INTEGER(iwp) ::  t_index        !<
123    INTEGER(iwp) ::  t_index_number !<
124   
125    LOGICAL  ::  reflect_x   !<
126    LOGICAL  ::  reflect_y   !<
127    LOGICAL  ::  reflect_z   !<
128
129    REAL(wp) ::  dt_particle !<
130    REAL(wp) ::  pos_x       !<
131    REAL(wp) ::  pos_x_old   !<
132    REAL(wp) ::  pos_y       !<
133    REAL(wp) ::  pos_y_old   !<
134    REAL(wp) ::  pos_z       !<
135    REAL(wp) ::  pos_z_old   !<
136    REAL(wp) ::  prt_x       !<
137    REAL(wp) ::  prt_y       !<
138    REAL(wp) ::  prt_z       !<
139    REAL(wp) ::  t(1:200)    !<
140    REAL(wp) ::  tmp_t       !<
141    REAL(wp) ::  xline       !<
142    REAL(wp) ::  yline       !<
143    REAL(wp) ::  zline       !<
144
145    IF ( range == 'bottom/top' )  THEN
146
147!
148!--    Apply boundary conditions to those particles that have crossed the top or
149!--    bottom boundary and delete those particles, which are older than allowed
150       DO  n = 1, number_of_particles
151
152          nn = particles(n)%tail_id
153
154!
155!--       Stop if particles have moved further than the length of one
156!--       PE subdomain (newly released particles have age = age_m!)
157          IF ( particles(n)%age /= particles(n)%age_m )  THEN
158             IF ( ABS(particles(n)%speed_x) >                                  &
159                  ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m)  .OR. &
160                  ABS(particles(n)%speed_y) >                                  &
161                  ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) )  THEN
162
163                  WRITE( message_string, * )  'particle too fast.  n = ',  n 
164                  CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
165             ENDIF
166          ENDIF
167
168          IF ( particles(n)%age > particle_maximum_age  .AND.  &
169               particles(n)%particle_mask )                              &
170          THEN
171             particles(n)%particle_mask  = .FALSE.
172             deleted_particles = deleted_particles + 1
173             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
174                tail_mask(nn) = .FALSE.
175                deleted_tails = deleted_tails + 1
176             ENDIF
177          ENDIF
178
179          IF ( particles(n)%z >= zu(nz)  .AND.  particles(n)%particle_mask )  THEN
180             IF ( ibc_par_t == 1 )  THEN
181!
182!--             Particle absorption
183                particles(n)%particle_mask  = .FALSE.
184                deleted_particles = deleted_particles + 1
185                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
186                   tail_mask(nn) = .FALSE.
187                   deleted_tails = deleted_tails + 1
188                ENDIF
189             ELSEIF ( ibc_par_t == 2 )  THEN
190!
191!--             Particle reflection
192                particles(n)%z       = 2.0_wp * zu(nz) - particles(n)%z
193                particles(n)%speed_z = -particles(n)%speed_z
194                IF ( use_sgs_for_particles  .AND. &
195                     particles(n)%rvar3 > 0.0_wp )  THEN
196                   particles(n)%rvar3 = -particles(n)%rvar3
197                ENDIF
198                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
199                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zu(nz) - &
200                                               particle_tail_coordinates(1,3,nn)
201                ENDIF
202             ENDIF
203          ENDIF
204         
205          IF ( particles(n)%z < zw(0)  .AND.  particles(n)%particle_mask )  THEN
206             IF ( ibc_par_b == 1 )  THEN
207!
208!--             Particle absorption
209                particles(n)%particle_mask  = .FALSE.
210                deleted_particles = deleted_particles + 1
211                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
212                   tail_mask(nn) = .FALSE.
213                   deleted_tails = deleted_tails + 1
214                ENDIF
215             ELSEIF ( ibc_par_b == 2 )  THEN
216!
217!--             Particle reflection
218                particles(n)%z       = 2.0_wp * zw(0) - particles(n)%z
219                particles(n)%speed_z = -particles(n)%speed_z
220                IF ( use_sgs_for_particles  .AND. &
221                     particles(n)%rvar3 < 0.0_wp )  THEN
222                   particles(n)%rvar3 = -particles(n)%rvar3
223                ENDIF
224                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
225                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zu(nz) - &
226                                               particle_tail_coordinates(1,3,nn)
227                ENDIF
228                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
229                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zw(0) - &
230                                               particle_tail_coordinates(1,3,nn)
231                ENDIF
232             ENDIF
233          ENDIF
234       ENDDO
235
236    ELSEIF ( range == 'walls' )  THEN
237
238       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
239
240       reflect_x = .FALSE.
241       reflect_y = .FALSE.
242       reflect_z = .FALSE.
243
244       DO  n = 1, number_of_particles
245
246          dt_particle = particles(n)%age - particles(n)%age_m
247
248          i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
249          j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
250          k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
251
252          prt_x = particles(n)%x
253          prt_y = particles(n)%y
254          prt_z = particles(n)%z
255
256!
257!--       If the particle position is below the surface, it has to be reflected
258          IF ( k2 <= nzb_s_inner(j2,i2)  .AND.  nzb_s_inner(j2,i2) /=0 )  THEN
259
260             pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
261             pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
262             pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
263             i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
264             j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
265             k1 = pos_z_old / dz + offset_ocean_nzt_m1
266
267!
268!--          Case 1
269             IF ( particles(n)%x > pos_x_old  .AND.  particles(n)%y > pos_y_old )&
270             THEN
271                t_index = 1
272
273                DO  i = i1, i2
274                   xline      = i * dx + 0.5_wp * dx
275                   t(t_index) = ( xline - pos_x_old ) / &
276                                ( particles(n)%x - pos_x_old )
277                   t_index    = t_index + 1
278                ENDDO
279
280                DO  j = j1, j2
281                   yline      = j * dy + 0.5_wp * dy
282                   t(t_index) = ( yline - pos_y_old ) / &
283                                ( particles(n)%y - pos_y_old )
284                   t_index    = t_index + 1
285                ENDDO
286
287                IF ( particles(n)%z < pos_z_old )  THEN
288                   DO  k = k1, k2 , -1
289                      zline      = k * dz
290                      t(t_index) = ( pos_z_old - zline ) / &
291                                   ( pos_z_old - particles(n)%z )
292                      t_index    = t_index + 1
293                   ENDDO
294                ENDIF
295
296                t_index_number = t_index - 1
297
298!
299!--             Sorting t
300                inc = 1
301                jr  = 1
302                DO WHILE ( inc <= t_index_number )
303                   inc = 3 * inc + 1
304                ENDDO
305
306                DO WHILE ( inc > 1 )
307                   inc = inc / 3
308                   DO  ir = inc+1, t_index_number
309                      tmp_t = t(ir)
310                      jr    = ir
311                      DO WHILE ( t(jr-inc) > tmp_t )
312                         t(jr) = t(jr-inc)
313                         jr    = jr - inc
314                         IF ( jr <= inc )  EXIT
315                      ENDDO
316                      t(jr) = tmp_t
317                   ENDDO
318                ENDDO
319
320         case1: DO  t_index = 1, t_index_number
321
322                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
323                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
324                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
325
326                   i3 = ( pos_x + 0.5_wp * dx ) * ddx   
327                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
328                   k3 = pos_z / dz + offset_ocean_nzt_m1
329
330                   i5 = pos_x * ddx
331                   j5 = pos_y * ddy
332                   k5 = pos_z / dz + offset_ocean_nzt_m1
333
334                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
335                        nzb_s_inner(j5,i3) /= 0 )  THEN
336
337                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
338                           k3 == nzb_s_inner(j5,i3) )  THEN
339                         reflect_z = .TRUE.
340                         EXIT case1
341                      ENDIF
342
343                   ENDIF
344
345                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
346                        nzb_s_inner(j3,i5) /= 0 )  THEN
347
348                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
349                           k3 == nzb_s_inner(j3,i5) )  THEN
350                         reflect_z = .TRUE.
351                         EXIT case1
352                      ENDIF
353
354                   ENDIF
355
356                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
357                        nzb_s_inner(j3,i3) /= 0 )  THEN
358
359                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
360                           k3 == nzb_s_inner(j3,i3) )  THEN
361                         reflect_z = .TRUE.
362                         EXIT case1
363                      ENDIF
364
365                      IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
366                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
367                         reflect_y = .TRUE.
368                         EXIT case1
369                      ENDIF
370
371                      IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
372                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
373                         reflect_x = .TRUE.
374                         EXIT case1
375                      ENDIF
376
377                   ENDIF
378
379                ENDDO case1
380
381!
382!--          Case 2
383             ELSEIF ( particles(n)%x > pos_x_old  .AND. &
384                      particles(n)%y < pos_y_old )  THEN
385
386                t_index = 1
387
388                DO  i = i1, i2
389                   xline      = i * dx + 0.5_wp * dx
390                   t(t_index) = ( xline - pos_x_old ) / &
391                                ( particles(n)%x - pos_x_old )
392                   t_index    = t_index + 1
393                ENDDO
394
395                DO  j = j1, j2, -1
396                   yline      = j * dy - 0.5_wp * dy
397                   t(t_index) = ( pos_y_old - yline ) / &
398                                ( pos_y_old - particles(n)%y )
399                   t_index    = t_index + 1
400                ENDDO
401
402                IF ( particles(n)%z < pos_z_old )  THEN
403                   DO  k = k1, k2 , -1
404                      zline      = k * dz
405                      t(t_index) = ( pos_z_old - zline ) / &
406                                   ( pos_z_old - particles(n)%z )
407                      t_index    = t_index + 1
408                   ENDDO
409                ENDIF
410                t_index_number = t_index-1
411
412!
413!--             Sorting t
414                inc = 1
415                jr  = 1
416                DO WHILE ( inc <= t_index_number )
417                   inc = 3 * inc + 1
418                ENDDO
419
420                DO WHILE ( inc > 1 )
421                   inc = inc / 3
422                   DO  ir = inc+1, t_index_number
423                      tmp_t = t(ir)
424                      jr    = ir
425                      DO WHILE ( t(jr-inc) > tmp_t )
426                         t(jr) = t(jr-inc)
427                         jr    = jr - inc
428                         IF ( jr <= inc )  EXIT
429                      ENDDO
430                      t(jr) = tmp_t
431                   ENDDO
432                ENDDO
433
434         case2: DO  t_index = 1, t_index_number
435
436                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
437                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
438                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
439
440                   i3 = ( pos_x + 0.5_wp * dx ) * ddx
441                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
442                   k3 = pos_z / dz + offset_ocean_nzt_m1
443
444                   i5 = pos_x * ddx
445                   j5 = pos_y * ddy
446                   k5 = pos_z / dz + offset_ocean_nzt_m1
447
448                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
449                        nzb_s_inner(j3,i5) /= 0 )  THEN
450
451                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
452                           k3 == nzb_s_inner(j3,i5) )  THEN
453                         reflect_z = .TRUE.
454                         EXIT case2
455                      ENDIF
456
457                   ENDIF
458
459                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
460                        nzb_s_inner(j3,i3) /= 0 )  THEN
461
462                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
463                           k3 == nzb_s_inner(j3,i3) )  THEN
464                         reflect_z = .TRUE.
465                         EXIT case2
466                      ENDIF
467
468                      IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
469                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
470                         reflect_x = .TRUE.
471                         EXIT case2
472                      ENDIF
473
474                   ENDIF
475
476                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
477                        nzb_s_inner(j5,i3) /= 0 )  THEN
478
479                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
480                           k3 == nzb_s_inner(j5,i3) )  THEN
481                         reflect_z = .TRUE.
482                         EXIT case2
483                      ENDIF
484
485                      IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
486                           pos_z < nzb_s_inner(j5,i3) * dz )  THEN
487                         reflect_y = .TRUE.
488                         EXIT case2
489                      ENDIF
490
491                   ENDIF
492
493                ENDDO case2
494
495!
496!--          Case 3
497             ELSEIF ( particles(n)%x < pos_x_old  .AND. &
498                      particles(n)%y > pos_y_old )  THEN
499
500                t_index = 1
501
502                DO  i = i1, i2, -1
503                   xline      = i * dx - 0.5_wp * dx
504                   t(t_index) = ( pos_x_old - xline ) / &
505                                ( pos_x_old - particles(n)%x )
506                   t_index    = t_index + 1
507                ENDDO
508
509                DO  j = j1, j2
510                   yline      = j * dy + 0.5_wp * dy
511                   t(t_index) = ( yline - pos_y_old ) / &
512                                ( particles(n)%y - pos_y_old )
513                   t_index    = t_index + 1
514                ENDDO
515
516                IF ( particles(n)%z < pos_z_old )  THEN
517                   DO  k = k1, k2 , -1
518                      zline      = k * dz
519                      t(t_index) = ( pos_z_old - zline ) / &
520                                   ( pos_z_old - particles(n)%z )
521                      t_index    = t_index + 1
522                   ENDDO
523                ENDIF
524                t_index_number = t_index - 1
525
526!
527!--             Sorting t
528                inc = 1
529                jr  = 1
530
531                DO WHILE ( inc <= t_index_number )
532                   inc = 3 * inc + 1
533                ENDDO
534
535                DO WHILE ( inc > 1 )
536                   inc = inc / 3
537                   DO  ir = inc+1, t_index_number
538                      tmp_t = t(ir)
539                      jr    = ir
540                      DO WHILE ( t(jr-inc) > tmp_t )
541                         t(jr) = t(jr-inc)
542                         jr    = jr - inc
543                         IF ( jr <= inc )  EXIT
544                      ENDDO
545                      t(jr) = tmp_t
546                   ENDDO
547                ENDDO
548
549         case3: DO  t_index = 1, t_index_number
550
551                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
552                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
553                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
554
555                   i3 = ( pos_x + 0.5_wp * dx ) * ddx
556                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
557                   k3 = pos_z / dz + offset_ocean_nzt_m1
558
559                   i5 = pos_x * ddx
560                   j5 = pos_y * ddy
561                   k5 = pos_z / dz + offset_ocean_nzt_m1
562
563                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
564                        nzb_s_inner(j5,i3) /= 0 )  THEN
565
566                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
567                           k3 == nzb_s_inner(j5,i3) )  THEN
568                         reflect_z = .TRUE.
569                         EXIT case3
570                      ENDIF
571
572                   ENDIF
573
574                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
575                        nzb_s_inner(j3,i3) /= 0 )  THEN
576
577                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
578                           k3 == nzb_s_inner(j3,i3) )  THEN
579                         reflect_z = .TRUE.
580                         EXIT case3
581                      ENDIF
582
583                      IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
584                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
585                         reflect_y = .TRUE.
586                         EXIT case3
587                      ENDIF
588
589                   ENDIF
590
591                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
592                        nzb_s_inner(j3,i5) /= 0 )  THEN
593
594                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
595                           k3 == nzb_s_inner(j3,i5) )  THEN
596                         reflect_z = .TRUE.
597                         EXIT case3
598                      ENDIF
599
600                      IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
601                           pos_z < nzb_s_inner(j3,i5) * dz )  THEN
602                         reflect_x = .TRUE.
603                         EXIT case3
604                      ENDIF
605
606                   ENDIF
607
608                ENDDO case3
609
610!
611!--          Case 4
612             ELSEIF ( particles(n)%x < pos_x_old  .AND. &
613                      particles(n)%y < pos_y_old )  THEN
614
615                t_index = 1
616
617                DO  i = i1, i2, -1
618                   xline      = i * dx - 0.5_wp * dx
619                   t(t_index) = ( pos_x_old - xline ) / &
620                                ( pos_x_old - particles(n)%x )
621                   t_index    = t_index + 1
622                ENDDO
623
624                DO  j = j1, j2, -1
625                   yline      = j * dy - 0.5_wp * dy
626                   t(t_index) = ( pos_y_old - yline ) / &
627                                ( pos_y_old - particles(n)%y )
628                   t_index    = t_index + 1
629                ENDDO
630
631                IF ( particles(n)%z < pos_z_old )  THEN
632                   DO  k = k1, k2 , -1
633                      zline      = k * dz
634                      t(t_index) = ( pos_z_old - zline ) / &
635                                   ( pos_z_old-particles(n)%z )
636                      t_index    = t_index + 1
637                   ENDDO
638                ENDIF
639                t_index_number = t_index-1
640
641!
642!--             Sorting t
643                inc = 1
644                jr  = 1
645
646                DO WHILE ( inc <= t_index_number )
647                   inc = 3 * inc + 1
648                ENDDO
649
650                DO WHILE ( inc > 1 )
651                   inc = inc / 3
652                   DO  ir = inc+1, t_index_number
653                      tmp_t = t(ir)
654                      jr = ir
655                      DO WHILE ( t(jr-inc) > tmp_t )
656                         t(jr) = t(jr-inc)
657                         jr    = jr - inc
658                         IF ( jr <= inc )  EXIT
659                      ENDDO
660                      t(jr) = tmp_t
661                   ENDDO
662                ENDDO
663
664         case4: DO  t_index = 1, t_index_number
665
666                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
667                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
668                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
669
670                   i3 = ( pos_x + 0.5_wp * dx ) * ddx   
671                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
672                   k3 = pos_z / dz + offset_ocean_nzt_m1
673
674                   i5 = pos_x * ddx
675                   j5 = pos_y * ddy
676                   k5 = pos_z / dz + offset_ocean_nzt_m1
677
678                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
679                        nzb_s_inner(j3,i3) /= 0 )  THEN
680
681                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
682                           k3 == nzb_s_inner(j3,i3) )  THEN
683                         reflect_z = .TRUE.
684                         EXIT case4
685                      ENDIF
686
687                   ENDIF
688
689                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
690                        nzb_s_inner(j3,i5) /= 0 )  THEN
691
692                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
693                           k3 == nzb_s_inner(j3,i5) )  THEN
694                         reflect_z = .TRUE.
695                         EXIT case4
696                      ENDIF
697
698                      IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
699                           nzb_s_inner(j3,i5) /=0  .AND.          &
700                           pos_z < nzb_s_inner(j3,i5) * dz )  THEN
701                         reflect_x = .TRUE.
702                         EXIT case4
703                      ENDIF
704
705                   ENDIF
706
707                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
708                        nzb_s_inner(j5,i3) /= 0 )  THEN
709
710                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
711                           k5 == nzb_s_inner(j5,i3) )  THEN
712                         reflect_z = .TRUE.
713                         EXIT case4
714                      ENDIF
715
716                      IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
717                           nzb_s_inner(j5,i3) /= 0  .AND.         &
718                           pos_z < nzb_s_inner(j5,i3) * dz )  THEN
719                         reflect_y = .TRUE.
720                         EXIT case4
721                      ENDIF
722
723                   ENDIF
724 
725                ENDDO case4
726
727             ENDIF
728
729          ENDIF      ! Check, if particle position is below the surface
730
731!
732!--       Do the respective reflection, in case that one of the above conditions
733!--       is found to be true
734          IF ( reflect_z )  THEN
735
736             particles(n)%z       = 2.0_wp * pos_z - prt_z
737             particles(n)%speed_z = - particles(n)%speed_z
738
739             IF ( use_sgs_for_particles )  THEN
740                particles(n)%rvar3 = - particles(n)%rvar3
741             ENDIF
742             reflect_z = .FALSE.
743
744          ELSEIF ( reflect_y )  THEN
745
746             particles(n)%y       = 2.0_wp * pos_y - prt_y
747             particles(n)%speed_y = - particles(n)%speed_y
748
749             IF ( use_sgs_for_particles )  THEN
750                particles(n)%rvar2 = - particles(n)%rvar2
751             ENDIF
752             reflect_y = .FALSE.
753
754          ELSEIF ( reflect_x )  THEN
755
756             particles(n)%x       = 2.0_wp * pos_x - prt_x
757             particles(n)%speed_x = - particles(n)%speed_x
758
759             IF ( use_sgs_for_particles )  THEN
760                particles(n)%rvar1 = - particles(n)%rvar1
761             ENDIF
762             reflect_x = .FALSE.
763
764          ENDIF       
765   
766       ENDDO
767
768       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
769
770    ENDIF
771
772 END SUBROUTINE lpm_boundary_conds
Note: See TracBrowser for help on using the repository browser.