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

Last change on this file since 1321 was 1321, checked in by raasch, 7 years ago

last commit documented

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