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

Last change on this file since 153 was 150, checked in by raasch, 16 years ago

particle advection allowed for ocean runs

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