Home

Context Navigation

source: palm/trunk/SOURCE/particle_boundary_conds.f90 @ 132

Last change on this file since 132 was 61, checked in by raasch, 18 years ago
further preliminary changes for revision 3.2
Property svn:keywords set to `Id`
File size: 17.9 KB

Rev	Line
[58]	1	SUBROUTINE particle_boundary_conds
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
	6	!
	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: particle_boundary_conds.f90 61 2007-03-12 05:42:06Z letzel $
	11	! Initial version (2007/03/09)
	12	!
	13	! Description:
	14	! ------------
	15	! Calculates the reflection of particles from vertical walls.
	16	! Routine developed by Jin Zhang (2006-2007)
[61]	17	! To do: Code structure for finding the t_index values and for checking the
	18	! reflection conditions is basically the same for all four cases, so it
	19	! should be possible to further simplify/shorten it.
[58]	20	!------------------------------------------------------------------------------!
	21
	22	USE control_parameters
	23	USE cpulog
	24	USE grid_variables
	25	USE indices
	26	USE interfaces
	27	USE particle_attributes
	28	USE pegrid
	29
	30	IMPLICIT NONE
	31
[60]	32	INTEGER :: i, inc, ir, i1, i2, i3, i5, j, jr, j1, j2, j3, j5, k, k1, k2, &
	33	k3, k5, n, t_index, t_index_number
[58]	34
[61]	35	LOGICAL :: reflect_x, reflect_y, reflect_z
	36
[60]	37	REAL :: dt_particle, pos_x, pos_x_old, pos_y, pos_y_old, pos_z, &
	38	pos_z_old, prt_x, prt_y, prt_z, tmp_t, xline, yline, zline
[58]	39
[60]	40	REAL :: t(1:200)
	41
[58]	42	CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' )
	43
	44
	45
[61]	46	reflect_x = .FALSE.
	47	reflect_y = .FALSE.
	48	reflect_z = .FALSE.
	49
[58]	50	DO n = 1, number_of_particles
	51
[60]	52	dt_particle = particles(n)%age - particles(n)%age_m
	53
[58]	54	i2 = ( particles(n)%x + 0.5 * dx ) * ddx
	55	j2 = ( particles(n)%y + 0.5 * dy ) * ddy
	56	k2 = particles(n)%z / dz + 1
	57
	58	prt_x = particles(n)%x
	59	prt_y = particles(n)%y
	60	prt_z = particles(n)%z
	61
[61]	62	!
	63	!-- If the particle position is below the surface, it has to be reflected
[58]	64	IF ( k2 <= nzb_s_inner(j2,i2) .AND. nzb_s_inner(j2,i2) /=0 ) THEN
	65
	66	pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
	67	pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
	68	pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
	69	i1 = ( pos_x_old + 0.5 * dx ) * ddx
	70	j1 = ( pos_y_old + 0.5 * dy ) * ddy
	71	k1 = pos_z_old / dz
	72
	73	!
	74	!-- Case 1
	75	IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old ) &
	76	THEN
	77	t_index = 1
	78
	79	DO i = i1, i2
	80	xline = i * dx + 0.5 * dx
	81	t(t_index) = ( xline - pos_x_old ) / &
	82	( particles(n)%x - pos_x_old )
	83	t_index = t_index + 1
	84	ENDDO
	85
	86	DO j = j1, j2
	87	yline = j * dy + 0.5 * dy
	88	t(t_index) = ( yline - pos_y_old ) / &
	89	( particles(n)%y - pos_y_old )
	90	t_index = t_index + 1
	91	ENDDO
	92
	93	IF ( particles(n)%z < pos_z_old ) THEN
	94	DO k = k1, k2 , -1
	95	zline = k * dz
	96	t(t_index) = ( pos_z_old - zline ) / &
	97	( pos_z_old - particles(n)%z )
	98	t_index = t_index + 1
	99	ENDDO
	100	ENDIF
	101
	102	t_index_number = t_index - 1
	103
	104	!
	105	!-- Sorting t
	106	inc = 1
	107	jr = 1
	108	DO WHILE ( inc <= t_index_number )
	109	inc = 3 * inc + 1
	110	ENDDO
	111
	112	DO WHILE ( inc > 1 )
	113	inc = inc / 3
	114	DO ir = inc+1, t_index_number
	115	tmp_t = t(ir)
	116	jr = ir
	117	DO WHILE ( t(jr-inc) > tmp_t )
	118	t(jr) = t(jr-inc)
	119	jr = jr - inc
	120	IF ( jr <= inc ) EXIT
	121	ENDDO
	122	t(jr) = tmp_t
	123	ENDDO
	124	ENDDO
	125
[61]	126	case1: DO t_index = 1, t_index_number
[58]	127
	128	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
	129	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
	130	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
	131
	132	i3 = ( pos_x + 0.5 * dx ) * ddx
	133	j3 = ( pos_y + 0.5 * dy ) * ddy
	134	k3 = pos_z / dz
	135
	136	i5 = pos_x * ddx
	137	j5 = pos_y * ddy
	138	k5 = pos_z / dz
	139
[61]	140	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
	141	nzb_s_inner(j5,i3) /= 0 ) THEN
	142
[58]	143	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
	144	k3 == nzb_s_inner(j5,i3) ) THEN
[61]	145	reflect_z = .TRUE.
	146	EXIT case1
	147	ENDIF
[58]	148
	149	ENDIF
	150
[61]	151	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
	152	nzb_s_inner(j3,i5) /= 0 ) THEN
	153
[58]	154	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
	155	k3 == nzb_s_inner(j3,i5) ) THEN
[61]	156	reflect_z = .TRUE.
	157	EXIT case1
	158	ENDIF
[58]	159
	160	ENDIF
	161
[61]	162	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
	163	nzb_s_inner(j3,i3) /= 0 ) THEN
	164
[58]	165	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
	166	k3 == nzb_s_inner(j3,i3) ) THEN
[61]	167	reflect_z = .TRUE.
	168	EXIT case1
[58]	169	ENDIF
	170
	171	IF ( pos_y == ( j3 * dy - 0.5 * dy ) .AND. &
	172	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
[61]	173	reflect_y = .TRUE.
	174	EXIT case1
[58]	175	ENDIF
	176
	177	IF ( pos_x == ( i3 * dx - 0.5 * dx ) .AND. &
	178	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
[61]	179	reflect_x = .TRUE.
	180	EXIT case1
	181	ENDIF
[58]	182
	183	ENDIF
	184
[61]	185	ENDDO case1
[58]	186
	187	!
	188	!-- Case 2
[61]	189	ELSEIF ( particles(n)%x > pos_x_old .AND. &
	190	particles(n)%y < pos_y_old ) THEN
	191
[58]	192	t_index = 1
	193
	194	DO i = i1, i2
	195	xline = i * dx + 0.5 * dx
	196	t(t_index) = ( xline - pos_x_old ) / &
	197	( particles(n)%x - pos_x_old )
	198	t_index = t_index + 1
	199	ENDDO
	200
	201	DO j = j1, j2, -1
	202	yline = j * dy - 0.5 * dy
	203	t(t_index) = ( pos_y_old - yline ) / &
	204	( pos_y_old - particles(n)%y )
	205	t_index = t_index + 1
	206	ENDDO
	207
	208	IF ( particles(n)%z < pos_z_old ) THEN
	209	DO k = k1, k2 , -1
	210	zline = k * dz
	211	t(t_index) = ( pos_z_old - zline ) / &
	212	( pos_z_old - particles(n)%z )
	213	t_index = t_index + 1
	214	ENDDO
	215	ENDIF
	216	t_index_number = t_index-1
	217
	218	!
	219	!-- Sorting t
	220	inc = 1
	221	jr = 1
	222	DO WHILE ( inc <= t_index_number )
	223	inc = 3 * inc + 1
	224	ENDDO
	225
	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
	237	ENDDO
	238	ENDDO
	239
[61]	240	case2: DO t_index = 1, t_index_number
[58]	241
	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 )
	245
	246	i3 = ( pos_x + 0.5 * dx ) * ddx
	247	j3 = ( pos_y + 0.5 * dy ) * ddy
	248	k3 = pos_z / dz
	249
	250	i5 = pos_x * ddx
	251	j5 = pos_y * ddy
	252	k5 = pos_z / dz
	253
[61]	254	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
	255	nzb_s_inner(j3,i5) /= 0 ) THEN
	256
[58]	257	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
	258	k3 == nzb_s_inner(j3,i5) ) THEN
[61]	259	reflect_z = .TRUE.
	260	EXIT case2
	261	ENDIF
[58]	262
	263	ENDIF
	264
[61]	265	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
	266	nzb_s_inner(j3,i3) /= 0 ) THEN
	267
[58]	268	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
	269	k3 == nzb_s_inner(j3,i3) ) THEN
[61]	270	reflect_z = .TRUE.
	271	EXIT case2
[58]	272	ENDIF
	273
	274	IF ( pos_x == ( i3 * dx - 0.5 * dx ) .AND. &
	275	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
[61]	276	reflect_x = .TRUE.
	277	EXIT case2
	278	ENDIF
[58]	279
	280	ENDIF
	281
[61]	282	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
	283	nzb_s_inner(j5,i3) /= 0 ) THEN
	284
[58]	285	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
	286	k3 == nzb_s_inner(j5,i3) ) THEN
[61]	287	reflect_z = .TRUE.
	288	EXIT case2
[58]	289	ENDIF
	290
	291	IF ( pos_y == ( j5 * dy + 0.5 * dy ) .AND. &
	292	pos_z < nzb_s_inner(j5,i3) * dz ) THEN
[61]	293	reflect_y = .TRUE.
	294	EXIT case2
	295	ENDIF
[58]	296
	297	ENDIF
	298
[61]	299	ENDDO case2
	300
[58]	301	!
	302	!-- Case 3
[61]	303	ELSEIF ( particles(n)%x < pos_x_old .AND. &
	304	particles(n)%y > pos_y_old ) THEN
	305
[58]	306	t_index = 1
	307
	308	DO i = i1, i2, -1
	309	xline = i * dx - 0.5 * dx
	310	t(t_index) = ( pos_x_old - xline ) / &
	311	( pos_x_old - particles(n)%x )
	312	t_index = t_index + 1
	313	ENDDO
	314
	315	DO j = j1, j2
	316	yline = j * dy + 0.5 * dy
	317	t(t_index) = ( yline - pos_y_old ) / &
	318	( particles(n)%y - pos_y_old )
	319	t_index = t_index + 1
	320	ENDDO
	321
	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
	332	!
	333	!-- Sorting t
	334	inc = 1
	335	jr = 1
	336
	337	DO WHILE ( inc <= t_index_number )
	338	inc = 3 * inc + 1
	339	ENDDO
	340
	341	DO WHILE ( inc > 1 )
	342	inc = inc / 3
	343	DO ir = inc+1, t_index_number
	344	tmp_t = t(ir)
	345	jr = ir
	346	DO WHILE ( t(jr-inc) > tmp_t )
	347	t(jr) = t(jr-inc)
	348	jr = jr - inc
	349	IF ( jr <= inc ) EXIT
	350	ENDDO
	351	t(jr) = tmp_t
	352	ENDDO
	353	ENDDO
	354
[61]	355	case3: DO t_index = 1, t_index_number
[58]	356
	357	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
	358	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
	359	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
	360
	361	i3 = ( pos_x + 0.5 * dx ) * ddx
	362	j3 = ( pos_y + 0.5 * dy ) * ddy
	363	k3 = pos_z / dz
	364
	365	i5 = pos_x * ddx
	366	j5 = pos_y * ddy
	367	k5 = pos_z / dz
	368
[61]	369	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
	370	nzb_s_inner(j5,i3) /= 0 ) THEN
	371
[58]	372	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
	373	k3 == nzb_s_inner(j5,i3) ) THEN
[61]	374	reflect_z = .TRUE.
	375	EXIT case3
	376	ENDIF
[58]	377
	378	ENDIF
	379
[61]	380	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
	381	nzb_s_inner(j3,i3) /= 0 ) THEN
	382
[58]	383	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
	384	k3 == nzb_s_inner(j3,i3) ) THEN
[61]	385	reflect_z = .TRUE.
	386	EXIT case3
[58]	387	ENDIF
	388
	389	IF ( pos_y == ( j3 * dy - 0.5 * dy ) .AND. &
	390	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
[61]	391	reflect_y = .TRUE.
	392	EXIT case3
	393	ENDIF
[58]	394
	395	ENDIF
	396
[61]	397	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
	398	nzb_s_inner(j3,i5) /= 0 ) THEN
	399
[58]	400	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
	401	k3 == nzb_s_inner(j3,i5) ) THEN
[61]	402	reflect_z = .TRUE.
	403	EXIT case3
[58]	404	ENDIF
	405
	406	IF ( pos_x == ( i5 * dx + 0.5 * dx ) .AND. &
	407	pos_z < nzb_s_inner(j3,i5) * dz ) THEN
[61]	408	reflect_x = .TRUE.
	409	EXIT case3
[58]	410	ENDIF
	411
[61]	412	ENDIF
	413
	414	ENDDO case3
	415
[58]	416	!
	417	!-- Case 4
[61]	418	ELSEIF ( particles(n)%x < pos_x_old .AND. &
	419	particles(n)%y < pos_y_old ) THEN
	420
[58]	421	t_index = 1
	422
	423	DO i = i1, i2, -1
	424	xline = i * dx - 0.5 * dx
	425	t(t_index) = ( pos_x_old - xline ) / &
	426	( pos_x_old - particles(n)%x )
	427	t_index = t_index + 1
	428	ENDDO
	429
	430	DO j = j1, j2, -1
	431	yline = j * dy - 0.5 * dy
	432	t(t_index) = ( pos_y_old - yline ) / &
	433	( pos_y_old - particles(n)%y )
	434	t_index = t_index + 1
	435	ENDDO
	436
	437	IF ( particles(n)%z < pos_z_old ) THEN
	438	DO k = k1, k2 , -1
	439	zline = k * dz
	440	t(t_index) = ( pos_z_old - zline ) / &
	441	( pos_z_old-particles(n)%z )
	442	t_index = t_index + 1
	443	ENDDO
	444	ENDIF
	445	t_index_number = t_index-1
	446
	447	!
	448	!-- Sorting t
	449	inc = 1
	450	jr = 1
	451
	452	DO WHILE ( inc <= t_index_number )
	453	inc = 3 * inc + 1
	454	ENDDO
	455
	456	DO WHILE ( inc > 1 )
	457	inc = inc / 3
	458	DO ir = inc+1, t_index_number
	459	tmp_t = t(ir)
	460	jr = ir
	461	DO WHILE ( t(jr-inc) > tmp_t )
	462	t(jr) = t(jr-inc)
	463	jr = jr - inc
	464	IF ( jr <= inc ) EXIT
	465	ENDDO
	466	t(jr) = tmp_t
	467	ENDDO
	468	ENDDO
	469
[61]	470	case4: DO t_index = 1, t_index_number
[58]	471
	472	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
	473	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
	474	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
	475
	476	i3 = ( pos_x + 0.5 * dx ) * ddx
	477	j3 = ( pos_y + 0.5 * dy ) * ddy
	478	k3 = pos_z / dz
	479
	480	i5 = pos_x * ddx
	481	j5 = pos_y * ddy
	482	k5 = pos_z / dz
	483
[61]	484	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
	485	nzb_s_inner(j3,i3) /= 0 ) THEN
	486
[58]	487	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
	488	k3 == nzb_s_inner(j3,i3) ) THEN
[61]	489	reflect_z = .TRUE.
	490	EXIT case4
	491	ENDIF
[58]	492
	493	ENDIF
	494
[61]	495	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
	496	nzb_s_inner(j3,i5) /= 0 ) THEN
	497
[58]	498	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
	499	k3 == nzb_s_inner(j3,i5) ) THEN
[61]	500	reflect_z = .TRUE.
	501	EXIT case4
[58]	502	ENDIF
	503
	504	IF ( pos_x == ( i5 * dx + 0.5 * dx ) .AND. &
	505	nzb_s_inner(j3,i5) /=0 .AND. &
	506	pos_z < nzb_s_inner(j3,i5) * dz ) THEN
[61]	507	reflect_x = .TRUE.
	508	EXIT case4
	509	ENDIF
[58]	510
	511	ENDIF
[61]	512
	513	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
	514	nzb_s_inner(j5,i3) /= 0 ) THEN
	515
[58]	516	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
	517	k5 == nzb_s_inner(j5,i3) ) THEN
[61]	518	reflect_z = .TRUE.
	519	EXIT case4
[58]	520	ENDIF
	521
	522	IF ( pos_y == ( j5 * dy + 0.5 * dy ) .AND. &
	523	nzb_s_inner(j5,i3) /= 0 .AND. &
	524	pos_z < nzb_s_inner(j5,i3) * dz ) THEN
[61]	525	reflect_y = .TRUE.
	526	EXIT case4
	527	ENDIF
[58]	528
[61]	529	ENDIF
	530
	531	ENDDO case4
	532
[58]	533	ENDIF
	534
[61]	535	ENDIF ! Check, if particle position is below the surface
	536
	537	!
	538	!-- Do the respective reflection, in case that one of the above conditions
	539	!-- is found to be true
	540	IF ( reflect_z ) THEN
	541
	542	particles(n)%z = 2.0 * pos_z - prt_z
	543	particles(n)%speed_z = - particles(n)%speed_z
	544
	545	IF ( use_sgs_for_particles ) THEN
	546	particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
	547	ENDIF
	548	reflect_z = .FALSE.
	549
	550	ELSEIF ( reflect_y ) THEN
	551
	552	particles(n)%y = 2.0 * pos_y - prt_y
	553	particles(n)%speed_y = - particles(n)%speed_y
	554
	555	IF ( use_sgs_for_particles ) THEN
	556	particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
	557	ENDIF
	558	reflect_y = .FALSE.
	559
	560	ELSEIF ( reflect_x ) THEN
	561
	562	particles(n)%x = 2.0 * pos_x - prt_x
	563	particles(n)%speed_x = - particles(n)%speed_x
	564
	565	IF ( use_sgs_for_particles ) THEN
	566	particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
	567	ENDIF
	568	reflect_x = .FALSE.
	569
	570	ENDIF
[58]	571
	572	ENDDO
	573
	574	CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' )
	575
[61]	576
[58]	577	END SUBROUTINE particle_boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |