Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |