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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

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