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

Last change on this file since 1359 was 1359, checked in by hoffmann, 7 years ago

new Lagrangian particle structure integrated

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