Home

Context Navigation

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

Last change on this file since 1319 was 1319, checked in by raasch, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 24.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |