Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |