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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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