Home

Context Navigation

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

Last change on this file since 1412 was 1360, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 17.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |