Home

Context Navigation

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

Last change on this file since 269 was 198, checked in by raasch, 16 years ago
file headers updated for the next release 3.5
Property svn:keywords set to `Id`
File size: 18.2 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 198 2008-09-17 08:55:28Z weinreis $
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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |