Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |