Home

Context Navigation

particle_boundary_conds.f90 @ 844

Last change on this file since 844 was 826, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 18.3 KB

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

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

Context Navigation

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

Download in other formats: