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

Last change on this file since 1905 was 1823, checked in by hoffmann, 9 years ago

last commit documented

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