Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |