Home

Context Navigation

lpm_boundary_conds.f90 @ 849

Last change on this file since 849 was 849, checked in by raasch, 12 years ago

Changed:

Original routine advec_particles split into several new subroutines and renamed
lpm.
init_particles renamed lpm_init
user_advec_particles renamed user_lpm_advec,
particle_boundary_conds renamed lpm_boundary_conds,
set_particle_attributes renamed lpm_set_attributes,
user_init_particles renamed user_lpm_init,
user_particle_attributes renamed user_lpm_set_attributes
(Makefile, lpm_droplet_collision, lpm_droplet_condensation, init_3d_model, modules, palm, read_var_list, time_integration, write_var_list, deleted: advec_particles, init_particles, particle_boundary_conds, set_particle_attributes, user_advec_particles, user_init_particles, user_particle_attributes, new: lpm, lpm_advec, lpm_boundary_conds, lpm_calc_liquid_water_content, lpm_data_output_particles, lpm_droplet_collision, lpm_drollet_condensation, lpm_exchange_horiz, lpm_extend_particle_array, lpm_extend_tails, lpm_extend_tail_array, lpm_init, lpm_init_sgs_tke, lpm_pack_arrays, lpm_read_restart_file, lpm_release_set, lpm_set_attributes, lpm_sort_arrays, lpm_write_exchange_statistics, lpm_write_restart_file, user_lpm_advec, user_lpm_init, user_lpm_set_attributes

Property svn:keywords set to Id

File size: 23.7 KB

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

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

Context Navigation

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

Download in other formats: