Home

Context Navigation

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

Last change on this file since 3709 was 3655, checked in by knoop, 6 years ago
Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update
Property svn:keywords set to `Id`
File size: 15.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |