Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |