Home

Context Navigation

source: palm/trunk/SOURCE/lpm_set_attributes.f90 @ 1827

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

Rev	Line
[1682]	1	!> @file lpm_set_attributes.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:
[264]	20	! -----------------
[1823]	21	!
[1360]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: lpm_set_attributes.f90 1823 2016-04-07 08:57:52Z maronga $
	26	!
[1823]	27	! 1822 2016-04-07 07:49:42Z hoffmann
	28	! 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
[264]	44	!
[1319]	45	! 1318 2014-03-17 13:35:16Z raasch
	46	! module interfaces removed
	47	!
[1037]	48	! 1036 2012-10-22 13:43:42Z raasch
	49	! code put under GPL (PALM 3.9)
	50	!
[850]	51	! 849 2012-03-15 10:35:09Z raasch
	52	! routine renamed: set_particle_attributes -> lpm_set_attributes
	53	!
[829]	54	! 828 2012-02-21 12:00:36Z raasch
	55	! particle feature color renamed class
	56	!
[392]	57	! 271 2009-03-26 00:47:14Z raasch
	58	! Initial version
[264]	59	!
	60	! Description:
	61	! ------------
[1682]	62	!> This routine sets certain particle attributes depending on the values that
	63	!> other PALM variables have at the current particle position.
[264]	64	!------------------------------------------------------------------------------!
[1682]	65	SUBROUTINE lpm_set_attributes
	66
[264]	67
[1320]	68	USE arrays_3d, &
	69	ONLY: pt, u, v, w, zu, zw
	70
	71	USE control_parameters, &
[1822]	72	ONLY: u_gtrans, v_gtrans, dz
[1320]	73
	74	USE cpulog, &
	75	ONLY: cpu_log, log_point_s
	76
	77	USE dvrp_variables, &
	78	ONLY: color_interval, dvrp_colortable_entries_prt, dvrpsize_interval, &
	79	particle_color, particle_dvrpsize
	80
	81	USE grid_variables, &
[1822]	82	ONLY: dx, dy
[1320]	83
	84	USE indices, &
	85	ONLY: ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt
	86
	87	USE kinds
	88
	89	USE particle_attributes, &
[1822]	90	ONLY: block_offset, grid_particles, number_of_particles, particles, &
	91	prt_count
[1320]	92
[264]	93	USE pegrid
	94
[1320]	95	USE statistics, &
	96	ONLY: sums, sums_l
	97
[264]	98	IMPLICIT NONE
	99
[1682]	100	INTEGER(iwp) :: i !<
	101	INTEGER(iwp) :: ip !<
	102	INTEGER(iwp) :: j !<
	103	INTEGER(iwp) :: jp !<
	104	INTEGER(iwp) :: k !<
	105	INTEGER(iwp) :: kp !<
	106	INTEGER(iwp) :: n !<
	107	INTEGER(iwp) :: nb !<
[264]	108
[1682]	109	INTEGER(iwp), DIMENSION(0:7) :: start_index !<
	110	INTEGER(iwp), DIMENSION(0:7) :: end_index !<
[1359]	111
[1682]	112	REAL(wp) :: aa !<
	113	REAL(wp) :: absuv !<
	114	REAL(wp) :: bb !<
	115	REAL(wp) :: cc !<
	116	REAL(wp) :: dd !<
	117	REAL(wp) :: gg !<
	118	REAL(wp) :: height !<
	119	REAL(wp) :: pt_int !<
	120	REAL(wp) :: pt_int_l !<
	121	REAL(wp) :: pt_int_u !<
	122	REAL(wp) :: u_int_l !<
	123	REAL(wp) :: u_int_u !<
	124	REAL(wp) :: v_int_l !<
	125	REAL(wp) :: v_int_u !<
	126	REAL(wp) :: w_int !<
	127	REAL(wp) :: w_int_l !<
	128	REAL(wp) :: w_int_u !<
	129	REAL(wp) :: x !<
	130	REAL(wp) :: y !<
[264]	131
[1682]	132	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_int !<
	133	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_int !<
	134	REAL(wp), DIMENSION(:), ALLOCATABLE :: xv !<
	135	REAL(wp), DIMENSION(:), ALLOCATABLE :: yv !<
	136	REAL(wp), DIMENSION(:), ALLOCATABLE :: zv !<
[1359]	137
[849]	138	CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'start' )
[264]	139
	140	!
	141	!-- Set particle color
	142	IF ( particle_color == 'absuv' ) THEN
	143
[266]	144	!
	145	!-- Set particle color depending on the absolute value of the horizontal
	146	!-- velocity
[1359]	147	DO ip = nxl, nxr
	148	DO jp = nys, nyn
	149	DO kp = nzb+1, nzt
[266]	150
[1359]	151	number_of_particles = prt_count(kp,jp,ip)
	152	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	153	IF ( number_of_particles <= 0 ) CYCLE
	154	start_index = grid_particles(kp,jp,ip)%start_index
	155	end_index = grid_particles(kp,jp,ip)%end_index
	156
	157	ALLOCATE( u_int(1:number_of_particles), &
	158	v_int(1:number_of_particles), &
	159	xv(1:number_of_particles), &
	160	yv(1:number_of_particles), &
	161	zv(1:number_of_particles) )
	162
	163	xv = particles(1:number_of_particles)%x
	164	yv = particles(1:number_of_particles)%y
	165	zv = particles(1:number_of_particles)%z
	166
	167	DO nb = 0,7
	168
	169	i = ip
	170	j = jp + block_offset(nb)%j_off
	171	k = kp + block_offset(nb)%k_off
	172
	173	DO n = start_index(nb), end_index(nb)
[266]	174	!
[1359]	175	!-- Interpolation of the velocity components in the xy-plane
	176	x = xv(n) + ( 0.5_wp - i ) * dx
	177	y = yv(n) - j * dy
	178	aa = x2 + y2
	179	bb = ( dx - x )2 + y2
	180	cc = x2 + ( dy - y )2
	181	dd = ( dx - x )2 + ( dy - y )2
	182	gg = aa + bb + cc + dd
[266]	183
[1359]	184	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * &
	185	u(k,j,i+1) + ( gg - cc ) * u(k,j+1,i) + &
	186	( gg - dd ) * u(k,j+1,i+1) &
	187	) / ( 3.0_wp * gg ) - u_gtrans
[266]	188
[1359]	189	IF ( k+1 == nzt+1 ) THEN
	190	u_int(n) = u_int_l
	191	ELSE
	192	u_int_u = ( ( gg - aa ) * u(k+1,j,i) + ( gg - bb ) * &
	193	u(k+1,j,i+1) + ( gg - cc ) * u(k+1,j+1,i) + &
	194	( gg - dd ) * u(k+1,j+1,i+1) &
	195	) / ( 3.0_wp * gg ) - u_gtrans
	196	u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz * &
	197	( u_int_u - u_int_l )
	198	ENDIF
	199
	200	ENDDO
	201
	202	i = ip + block_offset(nb)%i_off
	203	j = jp
	204	k = kp + block_offset(nb)%k_off
	205
	206	DO n = start_index(nb), end_index(nb)
[266]	207	!
[1359]	208	!-- Same procedure for interpolation of the v velocity-component
	209	x = xv(n) - i * dx
	210	y = yv(n) + ( 0.5_wp - j ) * dy
	211	aa = x2 + y2
	212	bb = ( dx - x )2 + y2
	213	cc = x2 + ( dy - y )2
	214	dd = ( dx - x )2 + ( dy - y )2
	215	gg = aa + bb + cc + dd
[266]	216
[1359]	217	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * &
	218	v(k,j,i+1) + ( gg - cc ) * v(k,j+1,i) + &
	219	( gg - dd ) * v(k,j+1,i+1) &
	220	) / ( 3.0_wp * gg ) - v_gtrans
[266]	221
[1359]	222	IF ( k+1 == nzt+1 ) THEN
	223	v_int(n) = v_int_l
	224	ELSE
	225	v_int_u = ( ( gg - aa ) * v(k+1,j,i) + ( gg - bb ) * &
	226	v(k+1,j,i+1) + ( gg - cc ) * v(k+1,j+1,i) + &
	227	( gg - dd ) * v(k+1,j+1,i+1) &
	228	) / ( 3.0_wp * gg ) - v_gtrans
	229	v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
	230	( v_int_u - v_int_l )
	231	ENDIF
[266]	232
[1359]	233	ENDDO
[266]	234
[1359]	235	ENDDO
	236
	237	DO n = 1, number_of_particles
	238
	239	absuv = SQRT( u_int(n)2 + v_int(n)2 )
	240
[266]	241	!
[1359]	242	!-- Limit values by the given interval and normalize to
	243	!-- interval [0,1]
	244	absuv = MIN( absuv, color_interval(2) )
	245	absuv = MAX( absuv, color_interval(1) )
[266]	246
[1359]	247	absuv = ( absuv - color_interval(1) ) / &
	248	( color_interval(2) - color_interval(1) )
[266]	249
	250	!
[1359]	251	!-- Number of available colors is defined in init_dvrp
	252	particles(n)%class = 1 + absuv * &
	253	( dvrp_colortable_entries_prt - 1 )
[266]	254
[1359]	255	ENDDO
	256
	257	DEALLOCATE( u_int, v_int, xv, yv, zv )
	258
	259	ENDDO
	260	ENDDO
[266]	261	ENDDO
	262
[264]	263	ELSEIF ( particle_color == 'pt*' ) THEN
	264	!
	265	!-- Set particle color depending on the resolved scale temperature
	266	!-- fluctuation.
	267	!-- First, calculate the horizontal average of the potential temperature
	268	!-- (This is also done in flow_statistics, but flow_statistics is called
	269	!-- after this routine.)
[1359]	270	sums_l(:,4,0) = 0.0_wp
[264]	271	DO i = nxl, nxr
	272	DO j = nys, nyn
	273	DO k = nzb, nzt+1
	274	sums_l(k,4,0) = sums_l(k,4,0) + pt(k,j,i)
	275	ENDDO
	276	ENDDO
	277	ENDDO
	278
	279	#if defined( __parallel )
[622]	280	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[264]	281	CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
	282	MPI_REAL, MPI_SUM, comm2d, ierr )
	283	#else
	284	sums(:,4) = sums_l(:,4,0)
	285	#endif
	286	sums(:,4) = sums(:,4) / ngp_2dh(0)
	287
[1359]	288	DO ip = nxl, nxr
	289	DO jp = nys, nyn
	290	DO kp = nzb+1, nzt
[264]	291
[1359]	292	number_of_particles = prt_count(kp,jp,ip)
	293	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	294	IF ( number_of_particles <= 0 ) CYCLE
	295	start_index = grid_particles(kp,jp,ip)%start_index
	296	end_index = grid_particles(kp,jp,ip)%end_index
[264]	297
[1359]	298	ALLOCATE( xv(1:number_of_particles), &
	299	yv(1:number_of_particles), &
	300	zv(1:number_of_particles) )
[264]	301
[1359]	302	xv = particles(1:number_of_particles)%x
	303	yv = particles(1:number_of_particles)%y
	304	zv = particles(1:number_of_particles)%z
[264]	305
[1359]	306	DO nb = 0,7
[264]	307
[1359]	308	i = ip + block_offset(nb)%i_off
	309	j = jp + block_offset(nb)%j_off
	310	k = kp + block_offset(nb)%k_off
	311
	312	DO n = start_index(nb), end_index(nb)
[264]	313	!
[1359]	314	!-- Interpolate temperature to the current particle position
	315	x = xv(n) - i * dx
	316	y = yv(n) - j * dy
	317	aa = x2 + y2
	318	bb = ( dx - x )2 + y2
	319	cc = x2 + ( dy - y )2
	320	dd = ( dx - x )2 + ( dy - y )2
	321	gg = aa + bb + cc + dd
[264]	322
[1359]	323	pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * &
	324	pt(k,j,i+1) + ( gg - cc ) * pt(k,j+1,i) + &
	325	( gg - dd ) * pt(k,j+1,i+1) &
	326	) / ( 3.0_wp * gg ) - sums(k,4)
[264]	327
[1359]	328	pt_int_u = ( ( gg - aa ) * pt(k+1,j,i) + ( gg - bb ) * &
	329	pt(k+1,j,i+1) + ( gg - cc ) * pt(k+1,j+1,i) + &
	330	( gg - dd ) * pt(k+1,j+1,i+1) &
	331	) / ( 3.0_wp * gg ) - sums(k,4)
	332
	333	pt_int = pt_int_l + ( zv(n) - zu(k) ) / dz * &
	334	( pt_int_u - pt_int_l )
	335
[264]	336	!
[1359]	337	!-- Limit values by the given interval and normalize to
	338	!-- interval [0,1]
	339	pt_int = MIN( pt_int, color_interval(2) )
	340	pt_int = MAX( pt_int, color_interval(1) )
[264]	341
[1359]	342	pt_int = ( pt_int - color_interval(1) ) / &
	343	( color_interval(2) - color_interval(1) )
	344
	345	!
	346	!-- Number of available colors is defined in init_dvrp
	347	particles(n)%class = 1 + pt_int * &
	348	( dvrp_colortable_entries_prt - 1 )
	349
	350	ENDDO
	351	ENDDO
	352
	353	DEALLOCATE( xv, yv, zv )
	354
	355	ENDDO
	356	ENDDO
[266]	357	ENDDO
[264]	358
	359	ELSEIF ( particle_color == 'z' ) THEN
[266]	360	!
	361	!-- Set particle color depending on the height above the bottom
	362	!-- boundary (z=0)
[1359]	363	DO ip = nxl, nxr
	364	DO jp = nys, nyn
	365	DO kp = nzb+1, nzt
[264]	366
[1359]	367	number_of_particles = prt_count(kp,jp,ip)
	368	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	369	IF ( number_of_particles <= 0 ) CYCLE
	370	DO n = 1, number_of_particles
	371
	372	height = particles(n)%z
[266]	373	!
[1359]	374	!-- Limit values by the given interval and normalize to
	375	!-- interval [0,1]
	376	height = MIN( height, color_interval(2) )
	377	height = MAX( height, color_interval(1) )
[266]	378
[1359]	379	height = ( height - color_interval(1) ) / &
	380	( color_interval(2) - color_interval(1) )
[266]	381
	382	!
[1359]	383	!-- Number of available colors is defined in init_dvrp
	384	particles(n)%class = 1 + height * &
	385	( dvrp_colortable_entries_prt - 1 )
[266]	386
[1359]	387	ENDDO
	388
	389	ENDDO
	390	ENDDO
[266]	391	ENDDO
	392
[264]	393	ENDIF
	394
	395	!
	396	!-- Set particle size for dvrp graphics
	397	IF ( particle_dvrpsize == 'absw' ) THEN
	398
[1359]	399	DO ip = nxl, nxr
	400	DO jp = nys, nyn
	401	DO kp = nzb+1, nzt
[264]	402
[1359]	403	number_of_particles = prt_count(kp,jp,ip)
	404	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	405	IF ( number_of_particles <= 0 ) CYCLE
	406	start_index = grid_particles(kp,jp,ip)%start_index
	407	end_index = grid_particles(kp,jp,ip)%end_index
[264]	408
[1359]	409	ALLOCATE( xv(1:number_of_particles), &
	410	yv(1:number_of_particles) )
[264]	411
[1359]	412	xv = particles(1:number_of_particles)%x
	413	yv = particles(1:number_of_particles)%y
	414	zv = particles(1:number_of_particles)%z
[264]	415
[1359]	416	DO nb = 0,7
	417
	418	i = ip + block_offset(nb)%i_off
	419	j = jp + block_offset(nb)%j_off
	420	k = kp-1
	421
	422	DO n = start_index(nb), end_index(nb)
[266]	423	!
[1359]	424	!-- Interpolate w-component to the current particle position
	425	x = xv(n) - i * dx
	426	y = yv(n) - j * dy
	427	aa = x2 + y2
	428	bb = ( dx - x )2 + y2
	429	cc = x2 + ( dy - y )2
	430	dd = ( dx - x )2 + ( dy - y )2
	431	gg = aa + bb + cc + dd
[264]	432
[1359]	433	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * &
	434	w(k,j,i+1) + ( gg - cc ) * w(k,j+1,i) + &
	435	( gg - dd ) * w(k,j+1,i+1) &
	436	) / ( 3.0_wp * gg )
[266]	437
[1359]	438	IF ( k+1 == nzt+1 ) THEN
	439	w_int = w_int_l
	440	ELSE
	441	w_int_u = ( ( gg - aa ) * w(k+1,j,i) + ( gg - bb ) * &
	442	w(k+1,j,i+1) + ( gg - cc ) * w(k+1,j+1,i) + &
	443	( gg - dd ) * w(k+1,j+1,i+1) &
	444	) / ( 3.0_wp * gg )
	445	w_int = w_int_l + ( zv(n) - zw(k) ) / dz * &
	446	( w_int_u - w_int_l )
	447	ENDIF
[266]	448
[1359]	449	!
	450	!-- Limit values by the given interval and normalize to
	451	!-- interval [0,1]
	452	w_int = ABS( w_int )
	453	w_int = MIN( w_int, dvrpsize_interval(2) )
	454	w_int = MAX( w_int, dvrpsize_interval(1) )
	455
	456	w_int = ( w_int - dvrpsize_interval(1) ) / &
	457	( dvrpsize_interval(2) - dvrpsize_interval(1) )
	458
	459	particles(n)%dvrp_psize = ( 0.25_wp + w_int * 0.6_wp ) * &
	460	dx
	461
	462	ENDDO
	463	ENDDO
	464
	465	DEALLOCATE( xv, yv, zv )
	466
	467	ENDDO
	468	ENDDO
[264]	469	ENDDO
	470
	471	ENDIF
	472
[849]	473	CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'stop' )
[264]	474
	475
[849]	476	END SUBROUTINE lpm_set_attributes

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |