Home

Context Navigation

source: palm/trunk/SOURCE/lpm_boundary_conds.f90 @ 1321

Last change on this file since 1321 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 26.8 KB

Line
1	SUBROUTINE lpm_boundary_conds( range )
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: lpm_boundary_conds.f90 1321 2014-03-20 09:40:40Z raasch $
26	!
27	! 1320 2014-03-20 08:40:49Z raasch
28	! ONLY-attribute added to USE-statements,
29	! kind-parameters added to all INTEGER and REAL declaration statements,
30	! kinds are defined in new module kinds,
31	! revision history before 2012 removed,
32	! comment fields (!:) to be used for variable explanations added to
33	! all variable declaration statements
34	!
35	! 1036 2012-10-22 13:43:42Z raasch
36	! code put under GPL (PALM 3.9)
37	!
38	! 849 2012-03-15 10:35:09Z raasch
39	! routine renamed lpm_boundary_conds, bottom and top boundary conditions
40	! included (former part of advec_particles)
41	!
42	! 824 2012-02-17 09:09:57Z raasch
43	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3
44	!
45	! Initial version (2007/03/09)
46	!
47	! Description:
48	! ------------
49	! Boundary conditions for the Lagrangian particles.
50	! The routine consists of two different parts. One handles the bottom (flat)
51	! and top boundary. In this part, also particles which exceeded their lifetime
52	! are deleted.
53	! The other part handles the reflection of particles from vertical walls.
54	! This part was developed by Jin Zhang during 2006-2007.
55	!
56	! To do: Code structure for finding the t_index values and for checking the
57	! ----- reflection conditions is basically the same for all four cases, so it
58	! should be possible to further simplify/shorten it.
59	!
60	! THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
61	! (see offset_ocean_*)
62	!------------------------------------------------------------------------------!
63
64	USE arrays_3d, &
65	ONLY: zu, zw
66
67	USE control_parameters, &
68	ONLY: dz, message_string, particle_maximum_age
69
70	USE cpulog, &
71	ONLY: cpu_log, log_point_s
72
73	USE grid_variables, &
74	ONLY: ddx, dx, ddy, dy
75
76	USE indices, &
77	ONLY: nxl, nxr, nyn, nys, nz, nzb_s_inner
78
79	USE kinds
80
81	USE particle_attributes, &
82	ONLY: deleted_particles, deleted_tails, ibc_par_b, ibc_par_t, &
83	number_of_particles, particles, particle_mask, &
84	particle_tail_coordinates, particle_type, offset_ocean_nzt_m1, &
85	tail_mask, use_particle_tails, use_sgs_for_particles
86
87	USE pegrid
88
89	IMPLICIT NONE
90
91	CHARACTER (LEN=*) :: range !:
92
93	INTEGER(iwp) :: i !:
94	INTEGER(iwp) :: inc !:
95	INTEGER(iwp) :: ir !:
96	INTEGER(iwp) :: i1 !:
97	INTEGER(iwp) :: i2 !:
98	INTEGER(iwp) :: i3 !:
99	INTEGER(iwp) :: i5 !:
100	INTEGER(iwp) :: j !:
101	INTEGER(iwp) :: jr !:
102	INTEGER(iwp) :: j1 !:
103	INTEGER(iwp) :: j2 !:
104	INTEGER(iwp) :: j3 !:
105	INTEGER(iwp) :: j5 !:
106	INTEGER(iwp) :: k !:
107	INTEGER(iwp) :: k1 !:
108	INTEGER(iwp) :: k2 !:
109	INTEGER(iwp) :: k3 !:
110	INTEGER(iwp) :: k5 !:
111	INTEGER(iwp) :: n !:
112	INTEGER(iwp) :: nn !:
113	INTEGER(iwp) :: t_index !:
114	INTEGER(iwp) :: t_index_number !:
115
116	LOGICAL :: reflect_x !:
117	LOGICAL :: reflect_y !:
118	LOGICAL :: reflect_z !:
119
120	REAL(wp) :: dt_particle !:
121	REAL(wp) :: pos_x !:
122	REAL(wp) :: pos_x_old !:
123	REAL(wp) :: pos_y !:
124	REAL(wp) :: pos_y_old !:
125	REAL(wp) :: pos_z !:
126	REAL(wp) :: pos_z_old !:
127	REAL(wp) :: prt_x !:
128	REAL(wp) :: prt_y !:
129	REAL(wp) :: prt_z !:
130	REAL(wp) :: t(1:200) !:
131	REAL(wp) :: tmp_t !:
132	REAL(wp) :: xline !:
133	REAL(wp) :: yline !:
134	REAL(wp) :: zline !:
135
136	IF ( range == 'bottom/top' ) THEN
137
138	!
139	!-- Apply boundary conditions to those particles that have crossed the top or
140	!-- bottom boundary and delete those particles, which are older than allowed
141	DO n = 1, number_of_particles
142
143	nn = particles(n)%tail_id
144
145	!
146	!-- Stop if particles have moved further than the length of one
147	!-- PE subdomain (newly released particles have age = age_m!)
148	IF ( particles(n)%age /= particles(n)%age_m ) THEN
149	IF ( ABS(particles(n)%speed_x) > &
150	((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. &
151	ABS(particles(n)%speed_y) > &
152	((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN
153
154	WRITE( message_string, * ) 'particle too fast. n = ', n
155	CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
156	ENDIF
157	ENDIF
158
159	IF ( particles(n)%age > particle_maximum_age .AND. &
160	particle_mask(n) ) &
161	THEN
162	particle_mask(n) = .FALSE.
163	deleted_particles = deleted_particles + 1
164	IF ( use_particle_tails .AND. nn /= 0 ) THEN
165	tail_mask(nn) = .FALSE.
166	deleted_tails = deleted_tails + 1
167	ENDIF
168	ENDIF
169
170	IF ( particles(n)%z >= zu(nz) .AND. particle_mask(n) ) THEN
171	IF ( ibc_par_t == 1 ) THEN
172	!
173	!-- Particle absorption
174	particle_mask(n) = .FALSE.
175	deleted_particles = deleted_particles + 1
176	IF ( use_particle_tails .AND. nn /= 0 ) THEN
177	tail_mask(nn) = .FALSE.
178	deleted_tails = deleted_tails + 1
179	ENDIF
180	ELSEIF ( ibc_par_t == 2 ) THEN
181	!
182	!-- Particle reflection
183	particles(n)%z = 2.0 * zu(nz) - particles(n)%z
184	particles(n)%speed_z = -particles(n)%speed_z
185	IF ( use_sgs_for_particles .AND. &
186	particles(n)%rvar3 > 0.0 ) THEN
187	particles(n)%rvar3 = -particles(n)%rvar3
188	ENDIF
189	IF ( use_particle_tails .AND. nn /= 0 ) THEN
190	particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
191	particle_tail_coordinates(1,3,nn)
192	ENDIF
193	ENDIF
194	ENDIF
195	IF ( particles(n)%z < zw(0) .AND. particle_mask(n) ) THEN
196	IF ( ibc_par_b == 1 ) THEN
197	!
198	!-- Particle absorption
199	particle_mask(n) = .FALSE.
200	deleted_particles = deleted_particles + 1
201	IF ( use_particle_tails .AND. nn /= 0 ) THEN
202	tail_mask(nn) = .FALSE.
203	deleted_tails = deleted_tails + 1
204	ENDIF
205	ELSEIF ( ibc_par_b == 2 ) THEN
206	!
207	!-- Particle reflection
208	particles(n)%z = 2.0 * zw(0) - particles(n)%z
209	particles(n)%speed_z = -particles(n)%speed_z
210	IF ( use_sgs_for_particles .AND. &
211	particles(n)%rvar3 < 0.0 ) THEN
212	particles(n)%rvar3 = -particles(n)%rvar3
213	ENDIF
214	IF ( use_particle_tails .AND. nn /= 0 ) THEN
215	particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
216	particle_tail_coordinates(1,3,nn)
217	ENDIF
218	IF ( use_particle_tails .AND. nn /= 0 ) THEN
219	particle_tail_coordinates(1,3,nn) = 2.0 * zw(0) - &
220	particle_tail_coordinates(1,3,nn)
221	ENDIF
222	ENDIF
223	ENDIF
224	ENDDO
225
226	ELSEIF ( range == 'walls' ) THEN
227
228	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
229
230	reflect_x = .FALSE.
231	reflect_y = .FALSE.
232	reflect_z = .FALSE.
233
234	DO n = 1, number_of_particles
235
236	dt_particle = particles(n)%age - particles(n)%age_m
237
238	i2 = ( particles(n)%x + 0.5 * dx ) * ddx
239	j2 = ( particles(n)%y + 0.5 * dy ) * ddy
240	k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
241
242	prt_x = particles(n)%x
243	prt_y = particles(n)%y
244	prt_z = particles(n)%z
245
246	!
247	!-- If the particle position is below the surface, it has to be reflected
248	IF ( k2 <= nzb_s_inner(j2,i2) .AND. nzb_s_inner(j2,i2) /=0 ) THEN
249
250	pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
251	pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
252	pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
253	i1 = ( pos_x_old + 0.5 * dx ) * ddx
254	j1 = ( pos_y_old + 0.5 * dy ) * ddy
255	k1 = pos_z_old / dz + offset_ocean_nzt_m1
256
257	!
258	!-- Case 1
259	IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old )&
260	THEN
261	t_index = 1
262
263	DO i = i1, i2
264	xline = i * dx + 0.5 * dx
265	t(t_index) = ( xline - pos_x_old ) / &
266	( particles(n)%x - pos_x_old )
267	t_index = t_index + 1
268	ENDDO
269
270	DO j = j1, j2
271	yline = j * dy + 0.5 * dy
272	t(t_index) = ( yline - pos_y_old ) / &
273	( particles(n)%y - pos_y_old )
274	t_index = t_index + 1
275	ENDDO
276
277	IF ( particles(n)%z < pos_z_old ) THEN
278	DO k = k1, k2 , -1
279	zline = k * dz
280	t(t_index) = ( pos_z_old - zline ) / &
281	( pos_z_old - particles(n)%z )
282	t_index = t_index + 1
283	ENDDO
284	ENDIF
285
286	t_index_number = t_index - 1
287
288	!
289	!-- Sorting t
290	inc = 1
291	jr = 1
292	DO WHILE ( inc <= t_index_number )
293	inc = 3 * inc + 1
294	ENDDO
295
296	DO WHILE ( inc > 1 )
297	inc = inc / 3
298	DO ir = inc+1, t_index_number
299	tmp_t = t(ir)
300	jr = ir
301	DO WHILE ( t(jr-inc) > tmp_t )
302	t(jr) = t(jr-inc)
303	jr = jr - inc
304	IF ( jr <= inc ) EXIT
305	ENDDO
306	t(jr) = tmp_t
307	ENDDO
308	ENDDO
309
310	case1: DO t_index = 1, t_index_number
311
312	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
313	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
314	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
315
316	i3 = ( pos_x + 0.5 * dx ) * ddx
317	j3 = ( pos_y + 0.5 * dy ) * ddy
318	k3 = pos_z / dz + offset_ocean_nzt_m1
319
320	i5 = pos_x * ddx
321	j5 = pos_y * ddy
322	k5 = pos_z / dz + offset_ocean_nzt_m1
323
324	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
325	nzb_s_inner(j5,i3) /= 0 ) THEN
326
327	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
328	k3 == nzb_s_inner(j5,i3) ) THEN
329	reflect_z = .TRUE.
330	EXIT case1
331	ENDIF
332
333	ENDIF
334
335	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
336	nzb_s_inner(j3,i5) /= 0 ) THEN
337
338	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
339	k3 == nzb_s_inner(j3,i5) ) THEN
340	reflect_z = .TRUE.
341	EXIT case1
342	ENDIF
343
344	ENDIF
345
346	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
347	nzb_s_inner(j3,i3) /= 0 ) THEN
348
349	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
350	k3 == nzb_s_inner(j3,i3) ) THEN
351	reflect_z = .TRUE.
352	EXIT case1
353	ENDIF
354
355	IF ( pos_y == ( j3 * dy - 0.5 * dy ) .AND. &
356	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
357	reflect_y = .TRUE.
358	EXIT case1
359	ENDIF
360
361	IF ( pos_x == ( i3 * dx - 0.5 * dx ) .AND. &
362	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
363	reflect_x = .TRUE.
364	EXIT case1
365	ENDIF
366
367	ENDIF
368
369	ENDDO case1
370
371	!
372	!-- Case 2
373	ELSEIF ( particles(n)%x > pos_x_old .AND. &
374	particles(n)%y < pos_y_old ) THEN
375
376	t_index = 1
377
378	DO i = i1, i2
379	xline = i * dx + 0.5 * dx
380	t(t_index) = ( xline - pos_x_old ) / &
381	( particles(n)%x - pos_x_old )
382	t_index = t_index + 1
383	ENDDO
384
385	DO j = j1, j2, -1
386	yline = j * dy - 0.5 * dy
387	t(t_index) = ( pos_y_old - yline ) / &
388	( pos_y_old - particles(n)%y )
389	t_index = t_index + 1
390	ENDDO
391
392	IF ( particles(n)%z < pos_z_old ) THEN
393	DO k = k1, k2 , -1
394	zline = k * dz
395	t(t_index) = ( pos_z_old - zline ) / &
396	( pos_z_old - particles(n)%z )
397	t_index = t_index + 1
398	ENDDO
399	ENDIF
400	t_index_number = t_index-1
401
402	!
403	!-- Sorting t
404	inc = 1
405	jr = 1
406	DO WHILE ( inc <= t_index_number )
407	inc = 3 * inc + 1
408	ENDDO
409
410	DO WHILE ( inc > 1 )
411	inc = inc / 3
412	DO ir = inc+1, t_index_number
413	tmp_t = t(ir)
414	jr = ir
415	DO WHILE ( t(jr-inc) > tmp_t )
416	t(jr) = t(jr-inc)
417	jr = jr - inc
418	IF ( jr <= inc ) EXIT
419	ENDDO
420	t(jr) = tmp_t
421	ENDDO
422	ENDDO
423
424	case2: DO t_index = 1, t_index_number
425
426	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
427	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
428	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
429
430	i3 = ( pos_x + 0.5 * dx ) * ddx
431	j3 = ( pos_y + 0.5 * dy ) * ddy
432	k3 = pos_z / dz + offset_ocean_nzt_m1
433
434	i5 = pos_x * ddx
435	j5 = pos_y * ddy
436	k5 = pos_z / dz + offset_ocean_nzt_m1
437
438	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
439	nzb_s_inner(j3,i5) /= 0 ) THEN
440
441	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
442	k3 == nzb_s_inner(j3,i5) ) THEN
443	reflect_z = .TRUE.
444	EXIT case2
445	ENDIF
446
447	ENDIF
448
449	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
450	nzb_s_inner(j3,i3) /= 0 ) THEN
451
452	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
453	k3 == nzb_s_inner(j3,i3) ) THEN
454	reflect_z = .TRUE.
455	EXIT case2
456	ENDIF
457
458	IF ( pos_x == ( i3 * dx - 0.5 * dx ) .AND. &
459	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
460	reflect_x = .TRUE.
461	EXIT case2
462	ENDIF
463
464	ENDIF
465
466	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
467	nzb_s_inner(j5,i3) /= 0 ) THEN
468
469	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
470	k3 == nzb_s_inner(j5,i3) ) THEN
471	reflect_z = .TRUE.
472	EXIT case2
473	ENDIF
474
475	IF ( pos_y == ( j5 * dy + 0.5 * dy ) .AND. &
476	pos_z < nzb_s_inner(j5,i3) * dz ) THEN
477	reflect_y = .TRUE.
478	EXIT case2
479	ENDIF
480
481	ENDIF
482
483	ENDDO case2
484
485	!
486	!-- Case 3
487	ELSEIF ( particles(n)%x < pos_x_old .AND. &
488	particles(n)%y > pos_y_old ) THEN
489
490	t_index = 1
491
492	DO i = i1, i2, -1
493	xline = i * dx - 0.5 * dx
494	t(t_index) = ( pos_x_old - xline ) / &
495	( pos_x_old - particles(n)%x )
496	t_index = t_index + 1
497	ENDDO
498
499	DO j = j1, j2
500	yline = j * dy + 0.5 * dy
501	t(t_index) = ( yline - pos_y_old ) / &
502	( particles(n)%y - pos_y_old )
503	t_index = t_index + 1
504	ENDDO
505
506	IF ( particles(n)%z < pos_z_old ) THEN
507	DO k = k1, k2 , -1
508	zline = k * dz
509	t(t_index) = ( pos_z_old - zline ) / &
510	( pos_z_old - particles(n)%z )
511	t_index = t_index + 1
512	ENDDO
513	ENDIF
514	t_index_number = t_index - 1
515
516	!
517	!-- Sorting t
518	inc = 1
519	jr = 1
520
521	DO WHILE ( inc <= t_index_number )
522	inc = 3 * inc + 1
523	ENDDO
524
525	DO WHILE ( inc > 1 )
526	inc = inc / 3
527	DO ir = inc+1, t_index_number
528	tmp_t = t(ir)
529	jr = ir
530	DO WHILE ( t(jr-inc) > tmp_t )
531	t(jr) = t(jr-inc)
532	jr = jr - inc
533	IF ( jr <= inc ) EXIT
534	ENDDO
535	t(jr) = tmp_t
536	ENDDO
537	ENDDO
538
539	case3: DO t_index = 1, t_index_number
540
541	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
542	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
543	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
544
545	i3 = ( pos_x + 0.5 * dx ) * ddx
546	j3 = ( pos_y + 0.5 * dy ) * ddy
547	k3 = pos_z / dz + offset_ocean_nzt_m1
548
549	i5 = pos_x * ddx
550	j5 = pos_y * ddy
551	k5 = pos_z / dz + offset_ocean_nzt_m1
552
553	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
554	nzb_s_inner(j5,i3) /= 0 ) THEN
555
556	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
557	k3 == nzb_s_inner(j5,i3) ) THEN
558	reflect_z = .TRUE.
559	EXIT case3
560	ENDIF
561
562	ENDIF
563
564	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
565	nzb_s_inner(j3,i3) /= 0 ) THEN
566
567	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
568	k3 == nzb_s_inner(j3,i3) ) THEN
569	reflect_z = .TRUE.
570	EXIT case3
571	ENDIF
572
573	IF ( pos_y == ( j3 * dy - 0.5 * dy ) .AND. &
574	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
575	reflect_y = .TRUE.
576	EXIT case3
577	ENDIF
578
579	ENDIF
580
581	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
582	nzb_s_inner(j3,i5) /= 0 ) THEN
583
584	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
585	k3 == nzb_s_inner(j3,i5) ) THEN
586	reflect_z = .TRUE.
587	EXIT case3
588	ENDIF
589
590	IF ( pos_x == ( i5 * dx + 0.5 * dx ) .AND. &
591	pos_z < nzb_s_inner(j3,i5) * dz ) THEN
592	reflect_x = .TRUE.
593	EXIT case3
594	ENDIF
595
596	ENDIF
597
598	ENDDO case3
599
600	!
601	!-- Case 4
602	ELSEIF ( particles(n)%x < pos_x_old .AND. &
603	particles(n)%y < pos_y_old ) THEN
604
605	t_index = 1
606
607	DO i = i1, i2, -1
608	xline = i * dx - 0.5 * dx
609	t(t_index) = ( pos_x_old - xline ) / &
610	( pos_x_old - particles(n)%x )
611	t_index = t_index + 1
612	ENDDO
613
614	DO j = j1, j2, -1
615	yline = j * dy - 0.5 * dy
616	t(t_index) = ( pos_y_old - yline ) / &
617	( pos_y_old - particles(n)%y )
618	t_index = t_index + 1
619	ENDDO
620
621	IF ( particles(n)%z < pos_z_old ) THEN
622	DO k = k1, k2 , -1
623	zline = k * dz
624	t(t_index) = ( pos_z_old - zline ) / &
625	( pos_z_old-particles(n)%z )
626	t_index = t_index + 1
627	ENDDO
628	ENDIF
629	t_index_number = t_index-1
630
631	!
632	!-- Sorting t
633	inc = 1
634	jr = 1
635
636	DO WHILE ( inc <= t_index_number )
637	inc = 3 * inc + 1
638	ENDDO
639
640	DO WHILE ( inc > 1 )
641	inc = inc / 3
642	DO ir = inc+1, t_index_number
643	tmp_t = t(ir)
644	jr = ir
645	DO WHILE ( t(jr-inc) > tmp_t )
646	t(jr) = t(jr-inc)
647	jr = jr - inc
648	IF ( jr <= inc ) EXIT
649	ENDDO
650	t(jr) = tmp_t
651	ENDDO
652	ENDDO
653
654	case4: DO t_index = 1, t_index_number
655
656	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
657	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
658	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
659
660	i3 = ( pos_x + 0.5 * dx ) * ddx
661	j3 = ( pos_y + 0.5 * dy ) * ddy
662	k3 = pos_z / dz + offset_ocean_nzt_m1
663
664	i5 = pos_x * ddx
665	j5 = pos_y * ddy
666	k5 = pos_z / dz + offset_ocean_nzt_m1
667
668	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
669	nzb_s_inner(j3,i3) /= 0 ) THEN
670
671	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
672	k3 == nzb_s_inner(j3,i3) ) THEN
673	reflect_z = .TRUE.
674	EXIT case4
675	ENDIF
676
677	ENDIF
678
679	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
680	nzb_s_inner(j3,i5) /= 0 ) THEN
681
682	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
683	k3 == nzb_s_inner(j3,i5) ) THEN
684	reflect_z = .TRUE.
685	EXIT case4
686	ENDIF
687
688	IF ( pos_x == ( i5 * dx + 0.5 * dx ) .AND. &
689	nzb_s_inner(j3,i5) /=0 .AND. &
690	pos_z < nzb_s_inner(j3,i5) * dz ) THEN
691	reflect_x = .TRUE.
692	EXIT case4
693	ENDIF
694
695	ENDIF
696
697	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
698	nzb_s_inner(j5,i3) /= 0 ) THEN
699
700	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
701	k5 == nzb_s_inner(j5,i3) ) THEN
702	reflect_z = .TRUE.
703	EXIT case4
704	ENDIF
705
706	IF ( pos_y == ( j5 * dy + 0.5 * dy ) .AND. &
707	nzb_s_inner(j5,i3) /= 0 .AND. &
708	pos_z < nzb_s_inner(j5,i3) * dz ) THEN
709	reflect_y = .TRUE.
710	EXIT case4
711	ENDIF
712
713	ENDIF
714
715	ENDDO case4
716
717	ENDIF
718
719	ENDIF ! Check, if particle position is below the surface
720
721	!
722	!-- Do the respective reflection, in case that one of the above conditions
723	!-- is found to be true
724	IF ( reflect_z ) THEN
725
726	particles(n)%z = 2.0 * pos_z - prt_z
727	particles(n)%speed_z = - particles(n)%speed_z
728
729	IF ( use_sgs_for_particles ) THEN
730	particles(n)%rvar3 = - particles(n)%rvar3
731	ENDIF
732	reflect_z = .FALSE.
733
734	ELSEIF ( reflect_y ) THEN
735
736	particles(n)%y = 2.0 * pos_y - prt_y
737	particles(n)%speed_y = - particles(n)%speed_y
738
739	IF ( use_sgs_for_particles ) THEN
740	particles(n)%rvar2 = - particles(n)%rvar2
741	ENDIF
742	reflect_y = .FALSE.
743
744	ELSEIF ( reflect_x ) THEN
745
746	particles(n)%x = 2.0 * pos_x - prt_x
747	particles(n)%speed_x = - particles(n)%speed_x
748
749	IF ( use_sgs_for_particles ) THEN
750	particles(n)%rvar1 = - particles(n)%rvar1
751	ENDIF
752	reflect_x = .FALSE.
753
754	ENDIF
755
756	ENDDO
757
758	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
759
760	ENDIF
761
762	END SUBROUTINE lpm_boundary_conds

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |