Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |