Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |