Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |