Home

Context Navigation

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

Last change on this file since 1929 was 1929, checked in by suehring, 8 years ago
several bugfixes in particle model and serial mode
Property svn:keywords set to `Id`
File size: 24.0 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	!
[1818]	16	! Copyright 1997-2016 Leibniz Universitaet Hannover
[1036]	17	!--------------------------------------------------------------------------------!
	18	!
[484]	19	! Current revisions:
[58]	20	! -----------------
[1929]	21	! Rewritten wall reflection
[1360]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: lpm_boundary_conds.f90 1929 2016-06-09 16:25:25Z suehring $
	26	!
[1823]	27	! 1822 2016-04-07 07:49:42Z hoffmann
	28	! Tails removed. Unused variables removed.
	29	!
[1683]	30	! 1682 2015-10-07 23:56:08Z knoop
	31	! Code annotations made doxygen readable
	32	!
[1360]	33	! 1359 2014-04-11 17:15:14Z hoffmann
	34	! New particle structure integrated.
	35	! Kind definition added to all floating point numbers.
	36	!
[1321]	37	! 1320 2014-03-20 08:40:49Z raasch
[1320]	38	! ONLY-attribute added to USE-statements,
	39	! kind-parameters added to all INTEGER and REAL declaration statements,
	40	! kinds are defined in new module kinds,
	41	! revision history before 2012 removed,
	42	! comment fields (!:) to be used for variable explanations added to
	43	! all variable declaration statements
[58]	44	!
[1037]	45	! 1036 2012-10-22 13:43:42Z raasch
	46	! code put under GPL (PALM 3.9)
	47	!
[850]	48	! 849 2012-03-15 10:35:09Z raasch
	49	! routine renamed lpm_boundary_conds, bottom and top boundary conditions
	50	! included (former part of advec_particles)
	51	!
[826]	52	! 824 2012-02-17 09:09:57Z raasch
	53	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3
	54	!
[58]	55	! Initial version (2007/03/09)
	56	!
	57	! Description:
	58	! ------------
[1682]	59	!> Boundary conditions for the Lagrangian particles.
	60	!> The routine consists of two different parts. One handles the bottom (flat)
	61	!> and top boundary. In this part, also particles which exceeded their lifetime
	62	!> are deleted.
	63	!> The other part handles the reflection of particles from vertical walls.
	64	!> This part was developed by Jin Zhang during 2006-2007.
	65	!>
	66	!> To do: Code structure for finding the t_index values and for checking the
	67	!> ----- reflection conditions is basically the same for all four cases, so it
	68	!> should be possible to further simplify/shorten it.
	69	!>
	70	!> THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
	71	!> (see offset_ocean_*)
[58]	72	!------------------------------------------------------------------------------!
[1682]	73	SUBROUTINE lpm_boundary_conds( range )
	74
[58]	75
[1320]	76	USE arrays_3d, &
	77	ONLY: zu, zw
[58]	78
[1320]	79	USE control_parameters, &
[1822]	80	ONLY: dz, message_string, particle_maximum_age
[58]	81
[1320]	82	USE cpulog, &
	83	ONLY: cpu_log, log_point_s
[849]	84
[1320]	85	USE grid_variables, &
	86	ONLY: ddx, dx, ddy, dy
[58]	87
[1320]	88	USE indices, &
	89	ONLY: nxl, nxr, nyn, nys, nz, nzb_s_inner
[61]	90
[1320]	91	USE kinds
[58]	92
[1320]	93	USE particle_attributes, &
[1822]	94	ONLY: deleted_particles, ibc_par_b, ibc_par_t, number_of_particles, &
	95	particles, particle_type, offset_ocean_nzt_m1, &
	96	use_sgs_for_particles
[60]	97
[1320]	98	USE pegrid
[58]	99
[1320]	100	IMPLICIT NONE
[58]	101
[1682]	102	CHARACTER (LEN=*) :: range !<
[1320]	103
[1929]	104	INTEGER(iwp) :: inc !< dummy for sorting algorithmus
	105	INTEGER(iwp) :: ir !< dummy for sorting algorithmus
	106	INTEGER(iwp) :: i1 !< grid index (x) of old particle position
	107	INTEGER(iwp) :: i2 !< grid index (x) of current particle position
	108	INTEGER(iwp) :: i3 !< grid index (x) of intermediate particle position
	109	INTEGER(iwp) :: jr !< dummy for sorting algorithmus
	110	INTEGER(iwp) :: j1 !< grid index (y) of old particle position
	111	INTEGER(iwp) :: j2 !< grid index (x) of current particle position
	112	INTEGER(iwp) :: j3 !< grid index (x) of intermediate particle position
	113	INTEGER(iwp) :: n !< particle number
	114	INTEGER(iwp) :: t_index !< running index for intermediate particle timesteps in reflection algorithmus
	115	INTEGER(iwp) :: t_index_number !< number of intermediate particle timesteps in reflection algorithmus
	116	INTEGER(iwp) :: tmp_x !< dummy for sorting algorithmus
	117	INTEGER(iwp) :: tmp_y !< dummy for sorting algorithmus
	118
	119	INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
	120	INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions
[1320]	121
[1929]	122	LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary
	123	LOGICAL :: cross_wall_y !< flag to check if particle reflection along y is necessary
	124	LOGICAL :: downwards !< flag to check if particle reflection along z is necessary (only if particle move downwards)
	125	LOGICAL :: reflect_x !< flag to check if particle is already reflected along x
	126	LOGICAL :: reflect_y !< flag to check if particle is already reflected along y
	127	LOGICAL :: reflect_z !< flag to check if particle is already reflected along z
	128	LOGICAL :: tmp_reach_x !< dummy for sorting algorithmus
	129	LOGICAL :: tmp_reach_y !< dummy for sorting algorithmus
	130	LOGICAL :: tmp_reach_z !< dummy for sorting algorithmus
	131	LOGICAL :: x_wall_reached !< flag to check if particle has already reached wall
	132	LOGICAL :: y_wall_reached !< flag to check if particle has already reached wall
[1320]	133
[1929]	134	LOGICAL, DIMENSION(0:10) :: reach_x !< flag to check if particle is at a yz-wall
	135	LOGICAL, DIMENSION(0:10) :: reach_y !< flag to check if particle is at a xz-wall
	136	LOGICAL, DIMENSION(0:10) :: reach_z !< flag to check if particle is at a xy-wall
[1320]	137
[1929]	138	REAL(wp) :: dt_particle !< particle timestep
	139	REAL(wp) :: dum !< dummy argument
	140	REAL(wp) :: eps = 1E-10_wp !< security number to check if particle has reached a wall
	141	REAL(wp) :: pos_x !< intermediate particle position (x)
	142	REAL(wp) :: pos_x_old !< particle position (x) at previous particle timestep
	143	REAL(wp) :: pos_y !< intermediate particle position (y)
	144	REAL(wp) :: pos_y_old !< particle position (y) at previous particle timestep
	145	REAL(wp) :: pos_z !< intermediate particle position (z)
	146	REAL(wp) :: pos_z_old !< particle position (z) at previous particle timestep
	147	REAL(wp) :: prt_x !< current particle position (x)
	148	REAL(wp) :: prt_y !< current particle position (y)
	149	REAL(wp) :: prt_z !< current particle position (z)
	150	REAL(wp) :: t_old !< previous reflection time
	151	REAL(wp) :: tmp_t !< dummy for sorting algorithmus
	152	REAL(wp) :: xwall !< location of wall in x
	153	REAL(wp) :: ywall !< location of wall in y
	154	REAL(wp) :: zwall1 !< location of wall in z (old grid box)
	155	REAL(wp) :: zwall2 !< location of wall in z (current grid box)
	156	REAL(wp) :: zwall3 !< location of wall in z (old y, current x)
	157	REAL(wp) :: zwall4 !< location of wall in z (current y, old x)
	158
	159	REAL(wp), DIMENSION(0:10) :: t !< reflection time
	160
	161
[849]	162	IF ( range == 'bottom/top' ) THEN
[58]	163
[849]	164	!
	165	!-- Apply boundary conditions to those particles that have crossed the top or
	166	!-- bottom boundary and delete those particles, which are older than allowed
	167	DO n = 1, number_of_particles
[61]	168
[849]	169	!
	170	!-- Stop if particles have moved further than the length of one
	171	!-- PE subdomain (newly released particles have age = age_m!)
	172	IF ( particles(n)%age /= particles(n)%age_m ) THEN
	173	IF ( ABS(particles(n)%speed_x) > &
	174	((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. &
	175	ABS(particles(n)%speed_y) > &
	176	((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN
[60]	177
[849]	178	WRITE( message_string, * ) 'particle too fast. n = ', n
	179	CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
	180	ENDIF
	181	ENDIF
[58]	182
[849]	183	IF ( particles(n)%age > particle_maximum_age .AND. &
[1359]	184	particles(n)%particle_mask ) &
[849]	185	THEN
[1359]	186	particles(n)%particle_mask = .FALSE.
[849]	187	deleted_particles = deleted_particles + 1
	188	ENDIF
[58]	189
[1359]	190	IF ( particles(n)%z >= zu(nz) .AND. particles(n)%particle_mask ) THEN
[849]	191	IF ( ibc_par_t == 1 ) THEN
[61]	192	!
[849]	193	!-- Particle absorption
[1359]	194	particles(n)%particle_mask = .FALSE.
[849]	195	deleted_particles = deleted_particles + 1
	196	ELSEIF ( ibc_par_t == 2 ) THEN
	197	!
	198	!-- Particle reflection
[1359]	199	particles(n)%z = 2.0_wp * zu(nz) - particles(n)%z
[849]	200	particles(n)%speed_z = -particles(n)%speed_z
	201	IF ( use_sgs_for_particles .AND. &
[1359]	202	particles(n)%rvar3 > 0.0_wp ) THEN
[849]	203	particles(n)%rvar3 = -particles(n)%rvar3
	204	ENDIF
	205	ENDIF
	206	ENDIF
[1359]	207
	208	IF ( particles(n)%z < zw(0) .AND. particles(n)%particle_mask ) THEN
[849]	209	IF ( ibc_par_b == 1 ) THEN
	210	!
	211	!-- Particle absorption
[1359]	212	particles(n)%particle_mask = .FALSE.
[849]	213	deleted_particles = deleted_particles + 1
	214	ELSEIF ( ibc_par_b == 2 ) THEN
	215	!
	216	!-- Particle reflection
[1359]	217	particles(n)%z = 2.0_wp * zw(0) - particles(n)%z
[849]	218	particles(n)%speed_z = -particles(n)%speed_z
	219	IF ( use_sgs_for_particles .AND. &
[1359]	220	particles(n)%rvar3 < 0.0_wp ) THEN
[849]	221	particles(n)%rvar3 = -particles(n)%rvar3
	222	ENDIF
	223	ENDIF
	224	ENDIF
	225	ENDDO
[58]	226
[849]	227	ELSEIF ( range == 'walls' ) THEN
[58]	228
[1929]	229
[849]	230	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
	231
	232	DO n = 1, number_of_particles
[1929]	233	!
	234	!-- Recalculate particle timestep
[849]	235	dt_particle = particles(n)%age - particles(n)%age_m
[1929]	236	!
	237	!-- Obtain x/y indices for current particle position
[1359]	238	i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
	239	j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
[1929]	240	!
	241	!-- Save current particle positions
[849]	242	prt_x = particles(n)%x
	243	prt_y = particles(n)%y
	244	prt_z = particles(n)%z
[58]	245	!
[1929]	246	!-- Recalculate old particle positions
	247	pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
	248	pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
	249	pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
[849]	250	!
[1929]	251	!-- Obtain x/y indices for old particle positions
	252	i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
	253	j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
[58]	254	!
[1929]	255	!-- Determine horizontal as well as vertical walls at which particle can
	256	!-- be potentially reflected.
	257	!-- Start with walls aligned in yz layer.
	258	!-- Wall to the right
	259	IF ( prt_x > pos_x_old ) THEN
	260	xwall = ( i1 + 0.5_wp ) * dx
	261	!
	262	!-- Wall to the left
	263	ELSE
	264	xwall = ( i1 - 0.5_wp ) * dx
	265	ENDIF
	266	!
	267	!-- Walls aligned in xz layer
	268	!-- Wall to the north
	269	IF ( prt_y > pos_y_old ) THEN
	270	ywall = ( j1 + 0.5_wp ) * dy
	271	!-- Wall to the south
	272	ELSE
	273	ywall = ( j1 - 0.5_wp ) * dy
	274	ENDIF
	275	!
	276	!-- Walls aligned in xy layer at which particle can be possiblly reflected
	277	zwall1 = zw(nzb_s_inner(j2,i2))
	278	zwall2 = zw(nzb_s_inner(j1,i1))
	279	zwall3 = zw(nzb_s_inner(j1,i2))
	280	zwall4 = zw(nzb_s_inner(j2,i1))
	281	!
	282	!-- Initialize flags to check if particle reflection is necessary
	283	downwards = .FALSE.
	284	cross_wall_x = .FALSE.
	285	cross_wall_y = .FALSE.
	286	!
	287	!-- Initialize flags to check if a wall is reached
	288	reach_x = .FALSE.
	289	reach_y = .FALSE.
	290	reach_z = .FALSE.
	291	!
	292	!-- Initialize flags to check if a particle was already reflected
	293	reflect_x = .FALSE.
	294	reflect_y = .FALSE.
	295	reflect_z = .FALSE.
	296	!
	297	!-- Initialize flags to check if a vertical wall is already crossed.
	298	!-- ( Required to obtain correct indices. )
	299	x_wall_reached = .FALSE.
	300	y_wall_reached = .FALSE.
	301	!
	302	!-- Initialize time array
	303	t = 0.0_wp
	304	!
	305	!-- Check if particle can reach any wall. This case, calculate the
	306	!-- fractional time needed to reach this wall. Store this fractional
	307	!-- timestep in array t. Moreover, store indices for these grid
	308	!-- boxes where the respective wall belongs to.
	309	!-- Start with x-direction.
	310	t_index = 1
	311	t(t_index) = ( xwall - pos_x_old ) &
	312	/ MERGE( MAX( prt_x - pos_x_old, 1E-30_wp ), &
	313	MIN( prt_x - pos_x_old, -1E-30_wp ), &
	314	prt_x > pos_x_old )
	315	x_ind(t_index) = i2
	316	y_ind(t_index) = j1
	317	reach_x(t_index) = .TRUE.
	318	reach_y(t_index) = .FALSE.
	319	reach_z(t_index) = .FALSE.
	320	!
	321	!-- Store these values only if particle really reaches any wall. t must
	322	!-- be in a interval between [0:1].
	323	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN
	324	t_index = t_index + 1
	325	cross_wall_x = .TRUE.
	326	ENDIF
	327	!
	328	!-- y-direction
	329	t(t_index) = ( ywall - pos_y_old ) &
	330	/ MERGE( MAX( prt_y - pos_y_old, 1E-30_wp ), &
	331	MIN( prt_y - pos_y_old, -1E-30_wp ), &
	332	prt_y > pos_y_old )
	333	x_ind(t_index) = i1
	334	y_ind(t_index) = j2
	335	reach_x(t_index) = .FALSE.
	336	reach_y(t_index) = .TRUE.
	337	reach_z(t_index) = .FALSE.
	338	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN
	339	t_index = t_index + 1
	340	cross_wall_y = .TRUE.
	341	ENDIF
	342	!
	343	!-- z-direction
	344	!-- At first, check if particle moves downwards. Only in this case a
	345	!-- particle can be reflected vertically.
	346	IF ( prt_z < pos_z_old ) THEN
[58]	347
[1929]	348	downwards = .TRUE.
	349	dum = 1.0_wp / MERGE( MAX( prt_z - pos_z_old, 1E-30_wp ), &
	350	MIN( prt_z - pos_z_old, -1E-30_wp ), &
	351	prt_z > pos_z_old )
[58]	352
[1929]	353	t(t_index) = ( zwall1 - pos_z_old ) * dum
	354	x_ind(t_index) = i2
	355	y_ind(t_index) = j2
	356	reach_x(t_index) = .FALSE.
	357	reach_y(t_index) = .FALSE.
	358	reach_z(t_index) = .TRUE.
	359	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) &
	360	t_index = t_index + 1
[58]	361
[1929]	362	reach_x(t_index) = .FALSE.
	363	reach_y(t_index) = .FALSE.
	364	reach_z(t_index) = .TRUE.
	365	t(t_index) = ( zwall2 - pos_z_old ) * dum
	366	x_ind(t_index) = i1
	367	y_ind(t_index) = j1
	368	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) &
	369	t_index = t_index + 1
[58]	370
[1929]	371	reach_x(t_index) = .FALSE.
	372	reach_y(t_index) = .FALSE.
	373	reach_z(t_index) = .TRUE.
	374	t(t_index) = ( zwall3 - pos_z_old ) * dum
	375	x_ind(t_index) = i2
	376	y_ind(t_index) = j1
	377	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) &
	378	t_index = t_index + 1
[58]	379
[1929]	380	reach_x(t_index) = .FALSE.
	381	reach_y(t_index) = .FALSE.
	382	reach_z(t_index) = .TRUE.
	383	t(t_index) = ( zwall4 - pos_z_old ) * dum
	384	x_ind(t_index) = i1
	385	y_ind(t_index) = j2
	386	IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) &
	387	t_index = t_index + 1
[58]	388
[1929]	389	END IF
	390	t_index_number = t_index - 1
[58]	391	!
[1929]	392	!-- Carry out reflection only if particle reaches any wall
	393	IF ( cross_wall_x .OR. cross_wall_y .OR. downwards ) THEN
[58]	394	!
[1929]	395	!-- Sort fractional timesteps in ascending order. Also sort the
	396	!-- corresponding indices and flag according to the time interval a
	397	!-- particle reaches the respective wall.
	398	inc = 1
	399	jr = 1
	400	DO WHILE ( inc <= t_index_number )
	401	inc = 3 * inc + 1
	402	ENDDO
[58]	403
[1929]	404	DO WHILE ( inc > 1 )
	405	inc = inc / 3
	406	DO ir = inc+1, t_index_number
	407	tmp_t = t(ir)
	408	tmp_x = x_ind(ir)
	409	tmp_y = y_ind(ir)
	410	tmp_reach_x = reach_x(ir)
	411	tmp_reach_y = reach_y(ir)
	412	tmp_reach_z = reach_z(ir)
	413	jr = ir
	414	DO WHILE ( t(jr-inc) > tmp_t )
	415	t(jr) = t(jr-inc)
	416	x_ind(jr) = x_ind(jr-inc)
	417	y_ind(jr) = y_ind(jr-inc)
	418	reach_x(jr) = reach_x(jr-inc)
	419	reach_y(jr) = reach_y(jr-inc)
	420	reach_z(jr) = reach_z(jr-inc)
	421	jr = jr - inc
	422	IF ( jr <= inc ) EXIT
[58]	423	ENDDO
[1929]	424	t(jr) = tmp_t
	425	x_ind(jr) = tmp_x
	426	y_ind(jr) = tmp_y
	427	reach_x(jr) = tmp_reach_x
	428	reach_y(jr) = tmp_reach_y
	429	reach_z(jr) = tmp_reach_z
[58]	430	ENDDO
[1929]	431	ENDDO
[58]	432	!
[1929]	433	!-- Initialize temporary particle positions
	434	pos_x = pos_x_old
	435	pos_y = pos_y_old
	436	pos_z = pos_z_old
	437	!
	438	!-- Loop over all times a particle possibly moves into a new grid box
	439	t_old = 0.0_wp
	440	DO t_index = 1, t_index_number
	441	!
	442	!-- Calculate intermediate particle position according to the
	443	!-- timesteps a particle reaches any wall.
	444	pos_x = pos_x + ( t(t_index) - t_old ) * dt_particle &
	445	* particles(n)%speed_x
	446	pos_y = pos_y + ( t(t_index) - t_old ) * dt_particle &
	447	* particles(n)%speed_y
	448	pos_z = pos_z + ( t(t_index) - t_old ) * dt_particle &
	449	* particles(n)%speed_z
	450	!
	451	!-- Obtain x/y grid indices for intermediate particle position from
	452	!-- sorted index array
	453	i3 = x_ind(t_index)
	454	j3 = y_ind(t_index)
	455	!
	456	!-- Check which wall is already reached
	457	IF ( .NOT. x_wall_reached ) x_wall_reached = reach_x(t_index)
	458	IF ( .NOT. y_wall_reached ) y_wall_reached = reach_y(t_index)
	459	!
	460	!-- Check if a particle needs to be reflected at any yz-wall. If
	461	!-- necessary, carry out reflection. Please note, a security
	462	!-- constant is required, as the particle position do not
	463	!-- necessarily exactly match the wall location due to rounding
	464	!-- errors.
	465	IF ( ABS( pos_x - xwall ) < eps .AND. &
	466	pos_z <= zw(nzb_s_inner(j3,i3)) .AND. &
	467	reach_x(t_index) .AND. &
	468	.NOT. reflect_x ) THEN
	469	!
	470	!-- Reflection in x-direction.
	471	!-- Ensure correct reflection by MIN/MAX functions, depending on
	472	!-- direction of particle transport.
	473	!-- Due to rounding errors pos_x do not exactly matches the wall
	474	!-- location, leading to erroneous reflection.
	475	pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ), &
	476	MAX( 2.0_wp * xwall - pos_x, xwall ), &
	477	particles(n)%x > xwall )
	478	!
	479	!-- Change sign of particle speed
	480	particles(n)%speed_x = - particles(n)%speed_x
	481	!
	482	!-- Change also sign of subgrid-scale particle speed
	483	particles(n)%rvar1 = - particles(n)%rvar1
	484	!
	485	!-- Set flag that reflection along x is already done
	486	reflect_x = .TRUE.
	487	!
	488	!-- As particle do not crosses any further yz-wall during
	489	!-- this timestep, set further x-indices to the current one.
	490	x_ind(t_index:t_index_number) = i1
	491	!
	492	!-- If particle already reached the wall but was not reflected,
	493	!-- set further x-indices to the new one.
	494	ELSEIF ( x_wall_reached .AND. .NOT. reflect_x ) THEN
	495	x_ind(t_index:t_index_number) = i2
	496	ENDIF
	497	!
	498	!-- Check if a particle needs to be reflected at any xz-wall. If
	499	!-- necessary, carry out reflection.
	500	IF ( ABS( pos_y - ywall ) < eps .AND. &
	501	pos_z <= zw(nzb_s_inner(j3,i3)) .AND. &
	502	reach_y(t_index) .AND. &
	503	.NOT. reflect_y ) THEN
[61]	504
[1929]	505	pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ), &
	506	MAX( 2.0_wp * ywall - pos_y, ywall ), &
	507	particles(n)%y > ywall )
[58]	508
[1929]	509	particles(n)%speed_y = - particles(n)%speed_y
	510	particles(n)%rvar2 = - particles(n)%rvar2
	511
	512	reflect_y = .TRUE.
	513	y_ind(t_index:t_index_number) = j1
[58]	514
[1929]	515	ELSEIF ( y_wall_reached .AND. .NOT. reflect_y ) THEN
	516	y_ind(t_index:t_index_number) = j2
[849]	517	ENDIF
[58]	518	!
[1929]	519	!-- Check if a particle needs to be reflected at any xy-wall. If
	520	!-- necessary, carry out reflection.
	521	IF ( downwards .AND. reach_z(t_index) .AND. &
	522	.NOT. reflect_z ) THEN
	523	IF ( pos_z - zw(nzb_s_inner(j3,i3)) < eps ) THEN
	524
	525	pos_z = MAX( 2.0_wp * zw(nzb_s_inner(j3,i3)) - pos_z, &
	526	zw(nzb_s_inner(j3,i3)) )
[58]	527
[1929]	528	particles(n)%speed_z = - particles(n)%speed_z
	529	particles(n)%rvar3 = - particles(n)%rvar3
[58]	530
[1929]	531	reflect_z = .TRUE.
[58]	532
[61]	533	ENDIF
[58]	534
[849]	535	ENDIF
[58]	536	!
[1929]	537	!-- Swap time
	538	t_old = t(t_index)
[58]	539
[1929]	540	ENDDO
[61]	541	!
[1929]	542	!-- If a particle was reflected, calculate final position from last
	543	!-- intermediate position.
	544	IF ( reflect_x .OR. reflect_y .OR. reflect_z ) THEN
[61]	545
[1929]	546	particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle &
	547	* particles(n)%speed_x
	548	particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle &
	549	* particles(n)%speed_y
	550	particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle &
	551	* particles(n)%speed_z
[61]	552
[849]	553	ENDIF
[61]	554
[1929]	555	ENDIF
[61]	556
[849]	557	ENDDO
[58]	558
[849]	559	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
[58]	560
[849]	561	ENDIF
[61]	562
[849]	563	END SUBROUTINE lpm_boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |