Home

Context Navigation

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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 27.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |