Home

Context Navigation

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

Last change on this file since 1834 was 1823, checked in by hoffmann, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 25.9 KB

Rev	Line
[1682]	1	!> @file lpm_boundary_conds.f90
[1036]	2	!--------------------------------------------------------------------------------!
	3	! This file is part of PALM.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the terms
	6	! of the GNU General Public License as published by the Free Software Foundation,
	7	! either version 3 of the License, or (at your option) any later version.
	8	!
	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	12	!
	13	! You should have received a copy of the GNU General Public License along with
	14	! PALM. If not, see <http://www.gnu.org/licenses/>.
	15	!
[1818]	16	! Copyright 1997-2016 Leibniz Universitaet Hannover
[1036]	17	!--------------------------------------------------------------------------------!
	18	!
[484]	19	! Current revisions:
[58]	20	! -----------------
[1823]	21	!
[1360]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: lpm_boundary_conds.f90 1823 2016-04-07 08:57:52Z raasch $
	26	!
[1823]	27	! 1822 2016-04-07 07:49:42Z hoffmann
	28	! Tails removed. Unused variables removed.
	29	!
[1683]	30	! 1682 2015-10-07 23:56:08Z knoop
	31	! Code annotations made doxygen readable
	32	!
[1360]	33	! 1359 2014-04-11 17:15:14Z hoffmann
	34	! New particle structure integrated.
	35	! Kind definition added to all floating point numbers.
	36	!
[1321]	37	! 1320 2014-03-20 08:40:49Z raasch
[1320]	38	! ONLY-attribute added to USE-statements,
	39	! kind-parameters added to all INTEGER and REAL declaration statements,
	40	! kinds are defined in new module kinds,
	41	! revision history before 2012 removed,
	42	! comment fields (!:) to be used for variable explanations added to
	43	! all variable declaration statements
[58]	44	!
[1037]	45	! 1036 2012-10-22 13:43:42Z raasch
	46	! code put under GPL (PALM 3.9)
	47	!
[850]	48	! 849 2012-03-15 10:35:09Z raasch
	49	! routine renamed lpm_boundary_conds, bottom and top boundary conditions
	50	! included (former part of advec_particles)
	51	!
[826]	52	! 824 2012-02-17 09:09:57Z raasch
	53	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3
	54	!
[58]	55	! Initial version (2007/03/09)
	56	!
	57	! Description:
	58	! ------------
[1682]	59	!> Boundary conditions for the Lagrangian particles.
	60	!> The routine consists of two different parts. One handles the bottom (flat)
	61	!> and top boundary. In this part, also particles which exceeded their lifetime
	62	!> are deleted.
	63	!> The other part handles the reflection of particles from vertical walls.
	64	!> This part was developed by Jin Zhang during 2006-2007.
	65	!>
	66	!> To do: Code structure for finding the t_index values and for checking the
	67	!> ----- reflection conditions is basically the same for all four cases, so it
	68	!> should be possible to further simplify/shorten it.
	69	!>
	70	!> THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
	71	!> (see offset_ocean_*)
[58]	72	!------------------------------------------------------------------------------!
[1682]	73	SUBROUTINE lpm_boundary_conds( range )
	74
[58]	75
[1320]	76	USE arrays_3d, &
	77	ONLY: zu, zw
[58]	78
[1320]	79	USE control_parameters, &
[1822]	80	ONLY: dz, message_string, particle_maximum_age
[58]	81
[1320]	82	USE cpulog, &
	83	ONLY: cpu_log, log_point_s
[849]	84
[1320]	85	USE grid_variables, &
	86	ONLY: ddx, dx, ddy, dy
[58]	87
[1320]	88	USE indices, &
	89	ONLY: nxl, nxr, nyn, nys, nz, nzb_s_inner
[61]	90
[1320]	91	USE kinds
[58]	92
[1320]	93	USE particle_attributes, &
[1822]	94	ONLY: deleted_particles, ibc_par_b, ibc_par_t, number_of_particles, &
	95	particles, particle_type, offset_ocean_nzt_m1, &
	96	use_sgs_for_particles
[60]	97
[1320]	98	USE pegrid
[58]	99
[1320]	100	IMPLICIT NONE
[58]	101
[1682]	102	CHARACTER (LEN=*) :: range !<
[1320]	103
[1682]	104	INTEGER(iwp) :: i !<
	105	INTEGER(iwp) :: inc !<
	106	INTEGER(iwp) :: ir !<
	107	INTEGER(iwp) :: i1 !<
	108	INTEGER(iwp) :: i2 !<
	109	INTEGER(iwp) :: i3 !<
	110	INTEGER(iwp) :: i5 !<
	111	INTEGER(iwp) :: j !<
	112	INTEGER(iwp) :: jr !<
	113	INTEGER(iwp) :: j1 !<
	114	INTEGER(iwp) :: j2 !<
	115	INTEGER(iwp) :: j3 !<
	116	INTEGER(iwp) :: j5 !<
	117	INTEGER(iwp) :: k !<
	118	INTEGER(iwp) :: k1 !<
	119	INTEGER(iwp) :: k2 !<
	120	INTEGER(iwp) :: k3 !<
	121	INTEGER(iwp) :: k5 !<
	122	INTEGER(iwp) :: n !<
	123	INTEGER(iwp) :: t_index !<
	124	INTEGER(iwp) :: t_index_number !<
[1320]	125
[1682]	126	LOGICAL :: reflect_x !<
	127	LOGICAL :: reflect_y !<
	128	LOGICAL :: reflect_z !<
[1320]	129
[1682]	130	REAL(wp) :: dt_particle !<
	131	REAL(wp) :: pos_x !<
	132	REAL(wp) :: pos_x_old !<
	133	REAL(wp) :: pos_y !<
	134	REAL(wp) :: pos_y_old !<
	135	REAL(wp) :: pos_z !<
	136	REAL(wp) :: pos_z_old !<
	137	REAL(wp) :: prt_x !<
	138	REAL(wp) :: prt_y !<
	139	REAL(wp) :: prt_z !<
	140	REAL(wp) :: t(1:200) !<
	141	REAL(wp) :: tmp_t !<
	142	REAL(wp) :: xline !<
	143	REAL(wp) :: yline !<
	144	REAL(wp) :: zline !<
[1320]	145
[849]	146	IF ( range == 'bottom/top' ) THEN
[58]	147
[849]	148	!
	149	!-- Apply boundary conditions to those particles that have crossed the top or
	150	!-- bottom boundary and delete those particles, which are older than allowed
	151	DO n = 1, number_of_particles
[61]	152
[849]	153	!
	154	!-- Stop if particles have moved further than the length of one
	155	!-- PE subdomain (newly released particles have age = age_m!)
	156	IF ( particles(n)%age /= particles(n)%age_m ) THEN
	157	IF ( ABS(particles(n)%speed_x) > &
	158	((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. &
	159	ABS(particles(n)%speed_y) > &
	160	((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN
[60]	161
[849]	162	WRITE( message_string, * ) 'particle too fast. n = ', n
	163	CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
	164	ENDIF
	165	ENDIF
[58]	166
[849]	167	IF ( particles(n)%age > particle_maximum_age .AND. &
[1359]	168	particles(n)%particle_mask ) &
[849]	169	THEN
[1359]	170	particles(n)%particle_mask = .FALSE.
[849]	171	deleted_particles = deleted_particles + 1
	172	ENDIF
[58]	173
[1359]	174	IF ( particles(n)%z >= zu(nz) .AND. particles(n)%particle_mask ) THEN
[849]	175	IF ( ibc_par_t == 1 ) THEN
[61]	176	!
[849]	177	!-- Particle absorption
[1359]	178	particles(n)%particle_mask = .FALSE.
[849]	179	deleted_particles = deleted_particles + 1
	180	ELSEIF ( ibc_par_t == 2 ) THEN
	181	!
	182	!-- Particle reflection
[1359]	183	particles(n)%z = 2.0_wp * zu(nz) - particles(n)%z
[849]	184	particles(n)%speed_z = -particles(n)%speed_z
	185	IF ( use_sgs_for_particles .AND. &
[1359]	186	particles(n)%rvar3 > 0.0_wp ) THEN
[849]	187	particles(n)%rvar3 = -particles(n)%rvar3
	188	ENDIF
	189	ENDIF
	190	ENDIF
[1359]	191
	192	IF ( particles(n)%z < zw(0) .AND. particles(n)%particle_mask ) THEN
[849]	193	IF ( ibc_par_b == 1 ) THEN
	194	!
	195	!-- Particle absorption
[1359]	196	particles(n)%particle_mask = .FALSE.
[849]	197	deleted_particles = deleted_particles + 1
	198	ELSEIF ( ibc_par_b == 2 ) THEN
	199	!
	200	!-- Particle reflection
[1359]	201	particles(n)%z = 2.0_wp * zw(0) - particles(n)%z
[849]	202	particles(n)%speed_z = -particles(n)%speed_z
	203	IF ( use_sgs_for_particles .AND. &
[1359]	204	particles(n)%rvar3 < 0.0_wp ) THEN
[849]	205	particles(n)%rvar3 = -particles(n)%rvar3
	206	ENDIF
	207	ENDIF
	208	ENDIF
	209	ENDDO
[58]	210
[849]	211	ELSEIF ( range == 'walls' ) THEN
[58]	212
[849]	213	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
	214
	215	reflect_x = .FALSE.
	216	reflect_y = .FALSE.
	217	reflect_z = .FALSE.
	218
	219	DO n = 1, number_of_particles
	220
	221	dt_particle = particles(n)%age - particles(n)%age_m
	222
[1359]	223	i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
	224	j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
[849]	225	k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
	226
	227	prt_x = particles(n)%x
	228	prt_y = particles(n)%y
	229	prt_z = particles(n)%z
	230
[58]	231	!
[849]	232	!-- If the particle position is below the surface, it has to be reflected
	233	IF ( k2 <= nzb_s_inner(j2,i2) .AND. nzb_s_inner(j2,i2) /=0 ) THEN
[58]	234
[849]	235	pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
	236	pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
	237	pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
[1359]	238	i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
	239	j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
[849]	240	k1 = pos_z_old / dz + offset_ocean_nzt_m1
[58]	241
[849]	242	!
	243	!-- Case 1
[1359]	244	IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old )&
[849]	245	THEN
	246	t_index = 1
[58]	247
[849]	248	DO i = i1, i2
[1359]	249	xline = i * dx + 0.5_wp * dx
[849]	250	t(t_index) = ( xline - pos_x_old ) / &
	251	( particles(n)%x - pos_x_old )
[58]	252	t_index = t_index + 1
	253	ENDDO
	254
[849]	255	DO j = j1, j2
[1359]	256	yline = j * dy + 0.5_wp * dy
[849]	257	t(t_index) = ( yline - pos_y_old ) / &
	258	( particles(n)%y - pos_y_old )
	259	t_index = t_index + 1
	260	ENDDO
[58]	261
[849]	262	IF ( particles(n)%z < pos_z_old ) THEN
	263	DO k = k1, k2 , -1
	264	zline = k * dz
	265	t(t_index) = ( pos_z_old - zline ) / &
	266	( pos_z_old - particles(n)%z )
	267	t_index = t_index + 1
	268	ENDDO
	269	ENDIF
	270
	271	t_index_number = t_index - 1
	272
[58]	273	!
[849]	274	!-- Sorting t
	275	inc = 1
	276	jr = 1
	277	DO WHILE ( inc <= t_index_number )
	278	inc = 3 * inc + 1
	279	ENDDO
[58]	280
[849]	281	DO WHILE ( inc > 1 )
	282	inc = inc / 3
	283	DO ir = inc+1, t_index_number
	284	tmp_t = t(ir)
	285	jr = ir
	286	DO WHILE ( t(jr-inc) > tmp_t )
	287	t(jr) = t(jr-inc)
	288	jr = jr - inc
	289	IF ( jr <= inc ) EXIT
	290	ENDDO
	291	t(jr) = tmp_t
[58]	292	ENDDO
	293	ENDDO
	294
[849]	295	case1: DO t_index = 1, t_index_number
[58]	296
[849]	297	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
	298	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
	299	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
[58]	300
[1359]	301	i3 = ( pos_x + 0.5_wp * dx ) * ddx
	302	j3 = ( pos_y + 0.5_wp * dy ) * ddy
[849]	303	k3 = pos_z / dz + offset_ocean_nzt_m1
[58]	304
[849]	305	i5 = pos_x * ddx
	306	j5 = pos_y * ddy
	307	k5 = pos_z / dz + offset_ocean_nzt_m1
[58]	308
[849]	309	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
	310	nzb_s_inner(j5,i3) /= 0 ) THEN
[61]	311
[849]	312	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
	313	k3 == nzb_s_inner(j5,i3) ) THEN
	314	reflect_z = .TRUE.
	315	EXIT case1
	316	ENDIF
	317
[61]	318	ENDIF
[58]	319
[849]	320	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
	321	nzb_s_inner(j3,i5) /= 0 ) THEN
[58]	322
[849]	323	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
	324	k3 == nzb_s_inner(j3,i5) ) THEN
	325	reflect_z = .TRUE.
	326	EXIT case1
	327	ENDIF
[61]	328
	329	ENDIF
[58]	330
[849]	331	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
	332	nzb_s_inner(j3,i3) /= 0 ) THEN
[58]	333
[849]	334	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
	335	k3 == nzb_s_inner(j3,i3) ) THEN
	336	reflect_z = .TRUE.
	337	EXIT case1
	338	ENDIF
[61]	339
[1359]	340	IF ( pos_y == ( j3 * dy - 0.5_wp * dy ) .AND. &
[849]	341	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
	342	reflect_y = .TRUE.
	343	EXIT case1
	344	ENDIF
[58]	345
[1359]	346	IF ( pos_x == ( i3 * dx - 0.5_wp * dx ) .AND. &
[849]	347	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
	348	reflect_x = .TRUE.
	349	EXIT case1
	350	ENDIF
[58]	351
[61]	352	ENDIF
[58]	353
[849]	354	ENDDO case1
[58]	355
	356	!
[849]	357	!-- Case 2
	358	ELSEIF ( particles(n)%x > pos_x_old .AND. &
	359	particles(n)%y < pos_y_old ) THEN
[61]	360
[849]	361	t_index = 1
[58]	362
[849]	363	DO i = i1, i2
[1359]	364	xline = i * dx + 0.5_wp * dx
[849]	365	t(t_index) = ( xline - pos_x_old ) / &
	366	( particles(n)%x - pos_x_old )
	367	t_index = t_index + 1
	368	ENDDO
[58]	369
[849]	370	DO j = j1, j2, -1
[1359]	371	yline = j * dy - 0.5_wp * dy
[849]	372	t(t_index) = ( pos_y_old - yline ) / &
	373	( pos_y_old - particles(n)%y )
[58]	374	t_index = t_index + 1
	375	ENDDO
	376
[849]	377	IF ( particles(n)%z < pos_z_old ) THEN
	378	DO k = k1, k2 , -1
	379	zline = k * dz
	380	t(t_index) = ( pos_z_old - zline ) / &
	381	( pos_z_old - particles(n)%z )
	382	t_index = t_index + 1
	383	ENDDO
	384	ENDIF
	385	t_index_number = t_index-1
	386
[58]	387	!
[849]	388	!-- Sorting t
	389	inc = 1
	390	jr = 1
	391	DO WHILE ( inc <= t_index_number )
	392	inc = 3 * inc + 1
	393	ENDDO
[58]	394
[849]	395	DO WHILE ( inc > 1 )
	396	inc = inc / 3
	397	DO ir = inc+1, t_index_number
	398	tmp_t = t(ir)
	399	jr = ir
	400	DO WHILE ( t(jr-inc) > tmp_t )
	401	t(jr) = t(jr-inc)
	402	jr = jr - inc
	403	IF ( jr <= inc ) EXIT
	404	ENDDO
	405	t(jr) = tmp_t
[58]	406	ENDDO
	407	ENDDO
	408
[849]	409	case2: DO t_index = 1, t_index_number
[58]	410
[849]	411	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
	412	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
	413	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
[58]	414
[1359]	415	i3 = ( pos_x + 0.5_wp * dx ) * ddx
	416	j3 = ( pos_y + 0.5_wp * dy ) * ddy
[849]	417	k3 = pos_z / dz + offset_ocean_nzt_m1
[58]	418
[849]	419	i5 = pos_x * ddx
	420	j5 = pos_y * ddy
	421	k5 = pos_z / dz + offset_ocean_nzt_m1
[58]	422
[849]	423	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
	424	nzb_s_inner(j3,i5) /= 0 ) THEN
[61]	425
[849]	426	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
	427	k3 == nzb_s_inner(j3,i5) ) THEN
	428	reflect_z = .TRUE.
	429	EXIT case2
	430	ENDIF
	431
[61]	432	ENDIF
[58]	433
[849]	434	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
	435	nzb_s_inner(j3,i3) /= 0 ) THEN
[58]	436
[849]	437	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
	438	k3 == nzb_s_inner(j3,i3) ) THEN
	439	reflect_z = .TRUE.
	440	EXIT case2
	441	ENDIF
[61]	442
[1359]	443	IF ( pos_x == ( i3 * dx - 0.5_wp * dx ) .AND. &
[849]	444	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
	445	reflect_x = .TRUE.
	446	EXIT case2
	447	ENDIF
[58]	448
[61]	449	ENDIF
[58]	450
[849]	451	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
	452	nzb_s_inner(j5,i3) /= 0 ) THEN
[58]	453
[849]	454	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
	455	k3 == nzb_s_inner(j5,i3) ) THEN
	456	reflect_z = .TRUE.
	457	EXIT case2
	458	ENDIF
[61]	459
[1359]	460	IF ( pos_y == ( j5 * dy + 0.5_wp * dy ) .AND. &
[849]	461	pos_z < nzb_s_inner(j5,i3) * dz ) THEN
	462	reflect_y = .TRUE.
	463	EXIT case2
	464	ENDIF
[58]	465
[61]	466	ENDIF
[58]	467
[849]	468	ENDDO case2
[58]	469
	470	!
[849]	471	!-- Case 3
	472	ELSEIF ( particles(n)%x < pos_x_old .AND. &
	473	particles(n)%y > pos_y_old ) THEN
[61]	474
[849]	475	t_index = 1
[58]	476
[849]	477	DO i = i1, i2, -1
[1359]	478	xline = i * dx - 0.5_wp * dx
[849]	479	t(t_index) = ( pos_x_old - xline ) / &
	480	( pos_x_old - particles(n)%x )
	481	t_index = t_index + 1
	482	ENDDO
[58]	483
[849]	484	DO j = j1, j2
[1359]	485	yline = j * dy + 0.5_wp * dy
[849]	486	t(t_index) = ( yline - pos_y_old ) / &
	487	( particles(n)%y - pos_y_old )
[58]	488	t_index = t_index + 1
	489	ENDDO
	490
[849]	491	IF ( particles(n)%z < pos_z_old ) THEN
	492	DO k = k1, k2 , -1
	493	zline = k * dz
	494	t(t_index) = ( pos_z_old - zline ) / &
	495	( pos_z_old - particles(n)%z )
	496	t_index = t_index + 1
	497	ENDDO
	498	ENDIF
	499	t_index_number = t_index - 1
	500
[58]	501	!
[849]	502	!-- Sorting t
	503	inc = 1
	504	jr = 1
[58]	505
[849]	506	DO WHILE ( inc <= t_index_number )
	507	inc = 3 * inc + 1
	508	ENDDO
[58]	509
[849]	510	DO WHILE ( inc > 1 )
	511	inc = inc / 3
	512	DO ir = inc+1, t_index_number
	513	tmp_t = t(ir)
	514	jr = ir
	515	DO WHILE ( t(jr-inc) > tmp_t )
	516	t(jr) = t(jr-inc)
	517	jr = jr - inc
	518	IF ( jr <= inc ) EXIT
	519	ENDDO
	520	t(jr) = tmp_t
[58]	521	ENDDO
	522	ENDDO
	523
[849]	524	case3: DO t_index = 1, t_index_number
[58]	525
[849]	526	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
	527	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
	528	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
[58]	529
[1359]	530	i3 = ( pos_x + 0.5_wp * dx ) * ddx
	531	j3 = ( pos_y + 0.5_wp * dy ) * ddy
[849]	532	k3 = pos_z / dz + offset_ocean_nzt_m1
[58]	533
[849]	534	i5 = pos_x * ddx
	535	j5 = pos_y * ddy
	536	k5 = pos_z / dz + offset_ocean_nzt_m1
[58]	537
[849]	538	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
	539	nzb_s_inner(j5,i3) /= 0 ) THEN
[61]	540
[849]	541	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
	542	k3 == nzb_s_inner(j5,i3) ) THEN
	543	reflect_z = .TRUE.
	544	EXIT case3
	545	ENDIF
	546
[61]	547	ENDIF
[58]	548
[849]	549	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
	550	nzb_s_inner(j3,i3) /= 0 ) THEN
[58]	551
[849]	552	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
	553	k3 == nzb_s_inner(j3,i3) ) THEN
	554	reflect_z = .TRUE.
	555	EXIT case3
	556	ENDIF
[61]	557
[1359]	558	IF ( pos_y == ( j3 * dy - 0.5_wp * dy ) .AND. &
[849]	559	pos_z < nzb_s_inner(j3,i3) * dz ) THEN
	560	reflect_y = .TRUE.
	561	EXIT case3
	562	ENDIF
[58]	563
[61]	564	ENDIF
[58]	565
[849]	566	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
	567	nzb_s_inner(j3,i5) /= 0 ) THEN
[58]	568
[849]	569	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
	570	k3 == nzb_s_inner(j3,i5) ) THEN
	571	reflect_z = .TRUE.
	572	EXIT case3
	573	ENDIF
[61]	574
[1359]	575	IF ( pos_x == ( i5 * dx + 0.5_wp * dx ) .AND. &
[849]	576	pos_z < nzb_s_inner(j3,i5) * dz ) THEN
	577	reflect_x = .TRUE.
	578	EXIT case3
	579	ENDIF
[58]	580
	581	ENDIF
	582
[849]	583	ENDDO case3
[61]	584
[58]	585	!
[849]	586	!-- Case 4
	587	ELSEIF ( particles(n)%x < pos_x_old .AND. &
	588	particles(n)%y < pos_y_old ) THEN
[61]	589
[849]	590	t_index = 1
[58]	591
[849]	592	DO i = i1, i2, -1
[1359]	593	xline = i * dx - 0.5_wp * dx
[849]	594	t(t_index) = ( pos_x_old - xline ) / &
	595	( pos_x_old - particles(n)%x )
	596	t_index = t_index + 1
	597	ENDDO
[58]	598
[849]	599	DO j = j1, j2, -1
[1359]	600	yline = j * dy - 0.5_wp * dy
[849]	601	t(t_index) = ( pos_y_old - yline ) / &
	602	( pos_y_old - particles(n)%y )
[58]	603	t_index = t_index + 1
	604	ENDDO
	605
[849]	606	IF ( particles(n)%z < pos_z_old ) THEN
	607	DO k = k1, k2 , -1
	608	zline = k * dz
	609	t(t_index) = ( pos_z_old - zline ) / &
	610	( pos_z_old-particles(n)%z )
	611	t_index = t_index + 1
	612	ENDDO
	613	ENDIF
	614	t_index_number = t_index-1
	615
[58]	616	!
[849]	617	!-- Sorting t
	618	inc = 1
	619	jr = 1
[58]	620
[849]	621	DO WHILE ( inc <= t_index_number )
	622	inc = 3 * inc + 1
	623	ENDDO
[58]	624
[849]	625	DO WHILE ( inc > 1 )
	626	inc = inc / 3
	627	DO ir = inc+1, t_index_number
	628	tmp_t = t(ir)
	629	jr = ir
	630	DO WHILE ( t(jr-inc) > tmp_t )
	631	t(jr) = t(jr-inc)
	632	jr = jr - inc
	633	IF ( jr <= inc ) EXIT
	634	ENDDO
	635	t(jr) = tmp_t
[58]	636	ENDDO
	637	ENDDO
	638
[849]	639	case4: DO t_index = 1, t_index_number
[58]	640
[849]	641	pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
	642	pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
	643	pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
[58]	644
[1359]	645	i3 = ( pos_x + 0.5_wp * dx ) * ddx
	646	j3 = ( pos_y + 0.5_wp * dy ) * ddy
[849]	647	k3 = pos_z / dz + offset_ocean_nzt_m1
[58]	648
[849]	649	i5 = pos_x * ddx
	650	j5 = pos_y * ddy
	651	k5 = pos_z / dz + offset_ocean_nzt_m1
[58]	652
[849]	653	IF ( k3 <= nzb_s_inner(j3,i3) .AND. &
	654	nzb_s_inner(j3,i3) /= 0 ) THEN
[61]	655
[849]	656	IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &
	657	k3 == nzb_s_inner(j3,i3) ) THEN
	658	reflect_z = .TRUE.
	659	EXIT case4
	660	ENDIF
	661
[61]	662	ENDIF
[58]	663
[849]	664	IF ( k5 <= nzb_s_inner(j3,i5) .AND. &
	665	nzb_s_inner(j3,i5) /= 0 ) THEN
[58]	666
[849]	667	IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &
	668	k3 == nzb_s_inner(j3,i5) ) THEN
	669	reflect_z = .TRUE.
	670	EXIT case4
	671	ENDIF
[61]	672
[1359]	673	IF ( pos_x == ( i5 * dx + 0.5_wp * dx ) .AND. &
[849]	674	nzb_s_inner(j3,i5) /=0 .AND. &
	675	pos_z < nzb_s_inner(j3,i5) * dz ) THEN
	676	reflect_x = .TRUE.
	677	EXIT case4
	678	ENDIF
[58]	679
[61]	680	ENDIF
[58]	681
[849]	682	IF ( k5 <= nzb_s_inner(j5,i3) .AND. &
	683	nzb_s_inner(j5,i3) /= 0 ) THEN
[61]	684
[849]	685	IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. &
	686	k5 == nzb_s_inner(j5,i3) ) THEN
	687	reflect_z = .TRUE.
	688	EXIT case4
	689	ENDIF
[61]	690
[1359]	691	IF ( pos_y == ( j5 * dy + 0.5_wp * dy ) .AND. &
[849]	692	nzb_s_inner(j5,i3) /= 0 .AND. &
	693	pos_z < nzb_s_inner(j5,i3) * dz ) THEN
	694	reflect_y = .TRUE.
	695	EXIT case4
	696	ENDIF
[58]	697
[61]	698	ENDIF
	699
[849]	700	ENDDO case4
[61]	701
[849]	702	ENDIF
[58]	703
[849]	704	ENDIF ! Check, if particle position is below the surface
[61]	705
	706	!
[849]	707	!-- Do the respective reflection, in case that one of the above conditions
	708	!-- is found to be true
	709	IF ( reflect_z ) THEN
[61]	710
[1359]	711	particles(n)%z = 2.0_wp * pos_z - prt_z
[849]	712	particles(n)%speed_z = - particles(n)%speed_z
[61]	713
[849]	714	IF ( use_sgs_for_particles ) THEN
	715	particles(n)%rvar3 = - particles(n)%rvar3
	716	ENDIF
	717	reflect_z = .FALSE.
[61]	718
[849]	719	ELSEIF ( reflect_y ) THEN
[61]	720
[1359]	721	particles(n)%y = 2.0_wp * pos_y - prt_y
[849]	722	particles(n)%speed_y = - particles(n)%speed_y
[61]	723
[849]	724	IF ( use_sgs_for_particles ) THEN
	725	particles(n)%rvar2 = - particles(n)%rvar2
	726	ENDIF
	727	reflect_y = .FALSE.
[61]	728
[849]	729	ELSEIF ( reflect_x ) THEN
[61]	730
[1359]	731	particles(n)%x = 2.0_wp * pos_x - prt_x
[849]	732	particles(n)%speed_x = - particles(n)%speed_x
[61]	733
[849]	734	IF ( use_sgs_for_particles ) THEN
	735	particles(n)%rvar1 = - particles(n)%rvar1
	736	ENDIF
	737	reflect_x = .FALSE.
[61]	738
[849]	739	ENDIF
[58]	740
[849]	741	ENDDO
[58]	742
[849]	743	CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
[58]	744
[849]	745	ENDIF
[61]	746
[849]	747	END SUBROUTINE lpm_boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |