Home

Context Navigation

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

Last change on this file since 1328 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 26.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |