Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |