Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |