Home

Context Navigation

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

Last change on this file since 2052 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 17.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |