Home

Context Navigation

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

Last change on this file since 161 was 150, checked in by raasch, 17 years ago
particle advection allowed for ocean runs
Property svn:keywords set to `Id`
File size: 18.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |