Home

Context Navigation

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

Last change on this file since 1629 was 1360, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 27.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |