Home

Context Navigation

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

Last change on this file since 2153 was 2101, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 24.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |