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

Last change on this file since 60 was 60, checked in by raasch, 15 years ago

preliminary update of further changes, running

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