Home

Context Navigation

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

Last change on this file since 60 was 60, checked in by raasch, 17 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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |