Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |