Home

Context Navigation

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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 27.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |