Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |