Home

Context Navigation

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

Last change on this file since 92 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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |