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

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

preliminary update of further changes, not running

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