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

Last change on this file since 550 was 484, checked in by raasch, 14 years ago

typo in file headers removed

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