Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |