Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |