Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |