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

Last change on this file since 143 was 61, checked in by raasch, 18 years ago

further preliminary changes for revision 3.2

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