Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |