source: palm/trunk/SOURCE/particle_boundary_conds.f90 @ 826

Last change on this file since 826 was 826, checked in by raasch, 10 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 18.3 KB
Line 
1 SUBROUTINE particle_boundary_conds
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: particle_boundary_conds.f90 826 2012-02-19 03:41:34Z raasch $
11!
12! 824 2012-02-17 09:09:57Z raasch
13! particle attributes speed_x|y|z_sgs renamed rvar1|2|3
14!
15! 150 2008-02-29 08:19:58Z raasch
16! Vertical index calculations adjusted for ocean runs.
17!
18! Initial version (2007/03/09)
19!
20! Description:
21! ------------
22! Calculates the reflection of particles from vertical walls.
23! Routine developed by Jin Zhang (2006-2007)
24! To do: Code structure for finding the t_index values and for checking the
25!        reflection conditions is basically the same for all four cases, so it
26!        should be possible to further simplify/shorten it.
27! THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR! (see offset_ocean_*)
28!------------------------------------------------------------------------------!
29
30    USE control_parameters
31    USE cpulog
32    USE grid_variables
33    USE indices
34    USE interfaces
35    USE particle_attributes
36    USE pegrid
37
38    IMPLICIT NONE
39
40    INTEGER ::  i, inc, ir, i1, i2, i3, i5, j, jr, j1, j2, j3, j5, k, k1, k2, &
41                k3, k5, n, t_index, t_index_number
42
43    LOGICAL ::  reflect_x, reflect_y, reflect_z
44
45    REAL ::  dt_particle, pos_x, pos_x_old, pos_y, pos_y_old, pos_z, &
46             pos_z_old, prt_x, prt_y, prt_z, tmp_t, xline, yline, zline
47
48    REAL ::  t(1:200)
49
50    CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' )
51
52
53
54    reflect_x = .FALSE.
55    reflect_y = .FALSE.
56    reflect_z = .FALSE.
57
58    DO  n = 1, number_of_particles
59
60       dt_particle = particles(n)%age - particles(n)%age_m
61
62       i2 = ( particles(n)%x + 0.5 * dx ) * ddx
63       j2 = ( particles(n)%y + 0.5 * dy ) * ddy
64       k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
65
66       prt_x = particles(n)%x
67       prt_y = particles(n)%y
68       prt_z = particles(n)%z
69
70!
71!--    If the particle position is below the surface, it has to be reflected
72       IF ( k2 <= nzb_s_inner(j2,i2)  .AND.  nzb_s_inner(j2,i2) /=0 )  THEN
73
74          pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
75          pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
76          pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
77          i1 = ( pos_x_old + 0.5 * dx ) * ddx
78          j1 = ( pos_y_old + 0.5 * dy ) * ddy
79          k1 = pos_z_old / dz + offset_ocean_nzt_m1
80
81!
82!--       Case 1
83          IF ( particles(n)%x > pos_x_old  .AND.  particles(n)%y > pos_y_old ) &
84          THEN
85             t_index = 1
86
87             DO  i = i1, i2
88                xline      = i * dx + 0.5 * dx
89                t(t_index) = ( xline - pos_x_old ) / &
90                             ( particles(n)%x - pos_x_old )
91                t_index    = t_index + 1
92             ENDDO
93
94             DO  j = j1, j2
95                yline      = j * dy + 0.5 * dy
96                t(t_index) = ( yline - pos_y_old ) / &
97                             ( particles(n)%y - pos_y_old )
98                t_index    = t_index + 1
99             ENDDO
100
101             IF ( particles(n)%z < pos_z_old )  THEN
102                DO  k = k1, k2 , -1
103                   zline      = k * dz
104                   t(t_index) = ( pos_z_old - zline ) / &
105                                ( pos_z_old - particles(n)%z )
106                   t_index    = t_index + 1
107                ENDDO
108             ENDIF
109
110             t_index_number = t_index - 1
111
112!
113!--          Sorting t
114             inc = 1
115             jr  = 1
116             DO WHILE ( inc <= t_index_number )
117                inc = 3 * inc + 1
118             ENDDO
119
120             DO WHILE ( inc > 1 )
121                inc = inc / 3
122                DO  ir = inc+1, t_index_number
123                   tmp_t = t(ir)
124                   jr    = ir
125                   DO WHILE ( t(jr-inc) > tmp_t )
126                      t(jr) = t(jr-inc)
127                      jr    = jr - inc
128                      IF ( jr <= inc )  EXIT
129                   ENDDO
130                   t(jr) = tmp_t
131                ENDDO
132             ENDDO
133
134      case1: DO  t_index = 1, t_index_number
135
136                pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
137                pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
138                pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
139
140                i3 = ( pos_x + 0.5 * dx ) * ddx   
141                j3 = ( pos_y + 0.5 * dy ) * ddy
142                k3 = pos_z / dz + offset_ocean_nzt_m1
143
144                i5 = pos_x * ddx
145                j5 = pos_y * ddy
146                k5 = pos_z / dz + offset_ocean_nzt_m1
147
148                IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
149                     nzb_s_inner(j5,i3) /= 0 )  THEN
150
151                   IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
152                        k3 == nzb_s_inner(j5,i3) )  THEN
153                      reflect_z = .TRUE.
154                      EXIT case1
155                   ENDIF
156
157                ENDIF
158
159                IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
160                     nzb_s_inner(j3,i5) /= 0 )  THEN
161
162                   IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
163                        k3 == nzb_s_inner(j3,i5) )  THEN
164                      reflect_z = .TRUE.
165                      EXIT case1
166                   ENDIF
167
168                ENDIF
169
170                IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
171                     nzb_s_inner(j3,i3) /= 0 )  THEN
172
173                   IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
174                        k3 == nzb_s_inner(j3,i3) )  THEN
175                      reflect_z = .TRUE.
176                      EXIT case1
177                   ENDIF
178
179                   IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
180                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
181                      reflect_y = .TRUE.
182                      EXIT case1
183                   ENDIF
184
185                   IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
186                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
187                      reflect_x = .TRUE.
188                      EXIT case1
189                   ENDIF
190
191                ENDIF
192
193             ENDDO case1
194
195!
196!--       Case 2
197          ELSEIF ( particles(n)%x > pos_x_old  .AND. &
198                   particles(n)%y < pos_y_old )  THEN
199
200             t_index = 1
201
202             DO  i = i1, i2
203                xline      = i * dx + 0.5 * dx
204                t(t_index) = ( xline - pos_x_old ) / &
205                             ( particles(n)%x - pos_x_old )
206                t_index    = t_index + 1
207             ENDDO
208
209             DO  j = j1, j2, -1
210                yline      = j * dy - 0.5 * dy
211                t(t_index) = ( pos_y_old - yline ) / &
212                             ( pos_y_old - particles(n)%y )
213                t_index    = t_index + 1
214             ENDDO
215
216             IF ( particles(n)%z < pos_z_old )  THEN
217                DO  k = k1, k2 , -1
218                   zline      = k * dz
219                   t(t_index) = ( pos_z_old - zline ) / &
220                                ( pos_z_old - particles(n)%z )
221                   t_index    = t_index + 1
222                ENDDO
223             ENDIF
224             t_index_number = t_index-1
225
226!
227!--          Sorting t
228             inc = 1
229             jr  = 1
230             DO WHILE ( inc <= t_index_number )
231                inc = 3 * inc + 1
232             ENDDO
233
234             DO WHILE ( inc > 1 )
235                inc = inc / 3
236                DO  ir = inc+1, t_index_number
237                   tmp_t = t(ir)
238                   jr    = ir
239                   DO WHILE ( t(jr-inc) > tmp_t )
240                      t(jr) = t(jr-inc)
241                      jr    = jr - inc
242                      IF ( jr <= inc )  EXIT
243                   ENDDO
244                   t(jr) = tmp_t
245                ENDDO
246             ENDDO
247
248      case2: DO  t_index = 1, t_index_number
249
250                pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
251                pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
252                pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
253
254                i3 = ( pos_x + 0.5 * dx ) * ddx
255                j3 = ( pos_y + 0.5 * dy ) * ddy
256                k3 = pos_z / dz + offset_ocean_nzt_m1
257
258                i5 = pos_x * ddx
259                j5 = pos_y * ddy
260                k5 = pos_z / dz + offset_ocean_nzt_m1
261
262                IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
263                     nzb_s_inner(j3,i5) /= 0 )  THEN
264
265                   IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
266                        k3 == nzb_s_inner(j3,i5) )  THEN
267                      reflect_z = .TRUE.
268                      EXIT case2
269                   ENDIF
270
271                ENDIF
272
273                IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
274                     nzb_s_inner(j3,i3) /= 0 )  THEN
275
276                   IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
277                        k3 == nzb_s_inner(j3,i3) )  THEN
278                      reflect_z = .TRUE.
279                      EXIT case2
280                   ENDIF
281
282                   IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
283                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
284                      reflect_x = .TRUE.
285                      EXIT case2
286                   ENDIF
287
288                ENDIF
289
290                IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
291                     nzb_s_inner(j5,i3) /= 0 )  THEN
292
293                   IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
294                        k3 == nzb_s_inner(j5,i3) )  THEN
295                      reflect_z = .TRUE.
296                      EXIT case2
297                   ENDIF
298
299                   IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
300                        pos_z < nzb_s_inner(j5,i3) * dz )  THEN
301                      reflect_y = .TRUE.
302                      EXIT case2
303                   ENDIF
304
305                ENDIF
306
307             ENDDO case2
308
309!
310!--       Case 3
311          ELSEIF ( particles(n)%x < pos_x_old  .AND. &
312                   particles(n)%y > pos_y_old )  THEN
313
314             t_index = 1
315
316             DO  i = i1, i2, -1
317                xline      = i * dx - 0.5 * dx
318                t(t_index) = ( pos_x_old - xline ) / &
319                             ( pos_x_old - particles(n)%x )
320                t_index    = t_index + 1
321             ENDDO
322
323             DO  j = j1, j2
324                yline      = j * dy + 0.5 * dy
325                t(t_index) = ( yline - pos_y_old ) / &
326                             ( particles(n)%y - pos_y_old )
327                t_index    = t_index + 1
328             ENDDO
329
330             IF ( particles(n)%z < pos_z_old )  THEN
331                DO  k = k1, k2 , -1
332                   zline      = k * dz
333                   t(t_index) = ( pos_z_old - zline ) / &
334                                ( pos_z_old - particles(n)%z )
335                   t_index    = t_index + 1
336                ENDDO
337             ENDIF
338             t_index_number = t_index - 1
339
340!
341!--          Sorting t
342             inc = 1
343             jr  = 1
344
345             DO WHILE ( inc <= t_index_number )
346                inc = 3 * inc + 1
347             ENDDO
348
349             DO WHILE ( inc > 1 )
350                inc = inc / 3
351                DO  ir = inc+1, t_index_number
352                   tmp_t = t(ir)
353                   jr    = ir
354                   DO WHILE ( t(jr-inc) > tmp_t )
355                      t(jr) = t(jr-inc)
356                      jr    = jr - inc
357                      IF ( jr <= inc )  EXIT
358                   ENDDO
359                   t(jr) = tmp_t
360                ENDDO
361             ENDDO
362
363      case3: DO  t_index = 1, t_index_number
364
365                pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
366                pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
367                pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
368
369                i3 = ( pos_x + 0.5 * dx ) * ddx
370                j3 = ( pos_y + 0.5 * dy ) * ddy
371                k3 = pos_z / dz + offset_ocean_nzt_m1
372
373                i5 = pos_x * ddx
374                j5 = pos_y * ddy
375                k5 = pos_z / dz + offset_ocean_nzt_m1
376
377                IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
378                     nzb_s_inner(j5,i3) /= 0 )  THEN
379
380                   IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
381                        k3 == nzb_s_inner(j5,i3) )  THEN
382                      reflect_z = .TRUE.
383                      EXIT case3
384                   ENDIF
385
386                ENDIF
387
388                IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
389                     nzb_s_inner(j3,i3) /= 0 )  THEN
390
391                   IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
392                        k3 == nzb_s_inner(j3,i3) )  THEN
393                      reflect_z = .TRUE.
394                      EXIT case3
395                   ENDIF
396
397                   IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
398                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
399                      reflect_y = .TRUE.
400                      EXIT case3
401                   ENDIF
402
403                ENDIF
404
405                IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
406                     nzb_s_inner(j3,i5) /= 0 )  THEN
407
408                   IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
409                        k3 == nzb_s_inner(j3,i5) )  THEN
410                      reflect_z = .TRUE.
411                      EXIT case3
412                   ENDIF
413
414                   IF ( pos_x == ( i5 * dx + 0.5 * dx )  .AND. &
415                        pos_z < nzb_s_inner(j3,i5) * dz )  THEN
416                      reflect_x = .TRUE.
417                      EXIT case3
418                   ENDIF
419
420                ENDIF
421
422             ENDDO case3
423
424!
425!--       Case 4
426          ELSEIF ( particles(n)%x < pos_x_old  .AND. &
427                   particles(n)%y < pos_y_old )  THEN
428
429             t_index = 1
430
431             DO  i = i1, i2, -1
432                xline      = i * dx - 0.5 * dx
433                t(t_index) = ( pos_x_old - xline ) / &
434                             ( pos_x_old - particles(n)%x )
435                t_index    = t_index + 1
436             ENDDO
437
438             DO  j = j1, j2, -1
439                yline      = j * dy - 0.5 * dy
440                t(t_index) = ( pos_y_old - yline ) / &
441                             ( pos_y_old - particles(n)%y )
442                t_index    = t_index + 1
443             ENDDO
444
445             IF ( particles(n)%z < pos_z_old )  THEN
446                DO  k = k1, k2 , -1
447                   zline      = k * dz
448                   t(t_index) = ( pos_z_old - zline ) / &
449                                ( pos_z_old-particles(n)%z )
450                   t_index    = t_index + 1
451                ENDDO
452             ENDIF
453             t_index_number = t_index-1
454
455!
456!--          Sorting t
457             inc = 1
458             jr  = 1
459
460             DO WHILE ( inc <= t_index_number )
461                inc = 3 * inc + 1
462             ENDDO
463
464             DO WHILE ( inc > 1 )
465                inc = inc / 3
466                DO  ir = inc+1, t_index_number
467                   tmp_t = t(ir)
468                   jr = ir
469                   DO WHILE ( t(jr-inc) > tmp_t )
470                      t(jr) = t(jr-inc)
471                      jr    = jr - inc
472                      IF ( jr <= inc )  EXIT
473                   ENDDO
474                   t(jr) = tmp_t
475                ENDDO
476             ENDDO
477
478      case4: DO  t_index = 1, t_index_number
479
480                pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
481                pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
482                pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
483
484                i3 = ( pos_x + 0.5 * dx ) * ddx   
485                j3 = ( pos_y + 0.5 * dy ) * ddy
486                k3 = pos_z / dz + offset_ocean_nzt_m1
487
488                i5 = pos_x * ddx
489                j5 = pos_y * ddy
490                k5 = pos_z / dz + offset_ocean_nzt_m1
491
492                IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
493                     nzb_s_inner(j3,i3) /= 0 )  THEN
494
495                   IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
496                        k3 == nzb_s_inner(j3,i3) )  THEN
497                      reflect_z = .TRUE.
498                      EXIT case4
499                   ENDIF
500
501                ENDIF
502
503                IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
504                     nzb_s_inner(j3,i5) /= 0 )  THEN
505
506                   IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
507                        k3 == nzb_s_inner(j3,i5) )  THEN
508                      reflect_z = .TRUE.
509                      EXIT case4
510                   ENDIF
511
512                   IF ( pos_x == ( i5 * dx + 0.5 * dx )  .AND. &
513                        nzb_s_inner(j3,i5) /=0  .AND.          &
514                        pos_z < nzb_s_inner(j3,i5) * dz )  THEN
515                      reflect_x = .TRUE.
516                      EXIT case4
517                   ENDIF
518
519                ENDIF
520
521                IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
522                     nzb_s_inner(j5,i3) /= 0 )  THEN
523
524                   IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
525                        k5 == nzb_s_inner(j5,i3) )  THEN
526                      reflect_z = .TRUE.
527                      EXIT case4
528                   ENDIF
529
530                   IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
531                        nzb_s_inner(j5,i3) /= 0  .AND.         &
532                        pos_z < nzb_s_inner(j5,i3) * dz )  THEN
533                      reflect_y = .TRUE.
534                      EXIT case4
535                   ENDIF
536
537                ENDIF
538 
539             ENDDO case4
540
541          ENDIF
542
543       ENDIF      ! Check, if particle position is below the surface
544
545!
546!--    Do the respective reflection, in case that one of the above conditions
547!--    is found to be true
548       IF ( reflect_z )  THEN
549
550          particles(n)%z       = 2.0 * pos_z - prt_z
551          particles(n)%speed_z = - particles(n)%speed_z
552
553          IF ( use_sgs_for_particles )  THEN
554             particles(n)%rvar3 = - particles(n)%rvar3
555          ENDIF
556          reflect_z = .FALSE.
557
558       ELSEIF ( reflect_y )  THEN
559
560          particles(n)%y       = 2.0 * pos_y - prt_y
561          particles(n)%speed_y = - particles(n)%speed_y
562
563          IF ( use_sgs_for_particles )  THEN
564             particles(n)%rvar2 = - particles(n)%rvar2
565          ENDIF
566          reflect_y = .FALSE.
567
568       ELSEIF ( reflect_x )  THEN
569
570          particles(n)%x       = 2.0 * pos_x - prt_x
571          particles(n)%speed_x = - particles(n)%speed_x
572
573          IF ( use_sgs_for_particles )  THEN
574             particles(n)%rvar1 = - particles(n)%rvar1
575          ENDIF
576          reflect_x = .FALSE.
577
578       ENDIF       
579   
580    ENDDO
581
582    CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' )
583
584
585 END SUBROUTINE particle_boundary_conds
Note: See TracBrowser for help on using the repository browser.