Home

Context Navigation

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

Last change on this file since 2707 was 2701, checked in by suehring, 7 years ago
changes from last commit documented
Property svn:keywords set to `Id`
File size: 27.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |