Home

Context Navigation

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |