Home

Context Navigation

source: palm/trunk/SOURCE/data_output_ptseries.f90 @ 1359

Last change on this file since 1359 was 1359, checked in by hoffmann, 10 years ago
new Lagrangian particle structure integrated
Property svn:keywords set to `Id`
File size: 16.3 KB

Rev	Line
[1]	1	SUBROUTINE data_output_ptseries
	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	!
[263]	20	! Current revisions:
[1]	21	! -----------------
[1359]	22	! New particle structure integrated.
[1329]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: data_output_ptseries.f90 1359 2014-04-11 17:15:14Z hoffmann $
	27	!
[1354]	28	! 1353 2014-04-08 15:21:23Z heinze
	29	! REAL constants provided with KIND-attribute
	30	!
[1329]	31	! 1327 2014-03-21 11:00:16Z raasch
	32	! -netcdf output queries
	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
[1]	41	!
[1319]	42	! 1318 2014-03-17 13:35:16Z raasch
	43	! barrier argument removed from cpu_log,
	44	! module interfaces removed
	45	!
[1037]	46	! 1036 2012-10-22 13:43:42Z raasch
	47	! code put under GPL (PALM 3.9)
	48	!
[826]	49	! 825 2012-02-19 03:03:44Z raasch
	50	! mean/minimum/maximum particle radius added as output quantity,
	51	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3
	52	!
[1]	53	! Revision 1.1 2006/08/04 14:24:18 raasch
	54	! Initial revision
	55	!
	56	!
	57	! Description:
	58	! ------------
	59	! Output of particle data timeseries in NetCDF format.
	60	!------------------------------------------------------------------------------!
	61
[1320]	62	USE cloud_parameters, &
	63	ONLY: curvature_solution_effects
	64
	65	USE control_parameters, &
[1327]	66	ONLY: dopts_time_count, time_since_reference_point
[1320]	67
	68	USE cpulog, &
	69	ONLY: cpu_log, log_point
	70
	71	USE indices, &
[1359]	72	ONLY: nxl, nxr, nys, nyn, nzb, nzt
[1320]	73
	74	USE kinds
	75
[1]	76	USE netcdf_control
[1320]	77
	78	USE particle_attributes, &
[1359]	79	ONLY: grid_particles, number_of_particles, number_of_particle_groups, &
	80	particles, prt_count
[1320]	81
[1]	82	USE pegrid
	83
	84	IMPLICIT NONE
	85
	86
[1320]	87	INTEGER(iwp) :: i !:
	88	INTEGER(iwp) :: inum !:
	89	INTEGER(iwp) :: j !:
[1359]	90	INTEGER(iwp) :: jg !:
	91	INTEGER(iwp) :: k !:
[1320]	92	INTEGER(iwp) :: n !:
[1]	93
[1320]	94	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pts_value !:
	95	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pts_value_l !:
[1]	96
	97
	98	CALL cpu_log( log_point(36), 'data_output_ptseries', 'start' )
	99
[1327]	100	IF ( myid == 0 ) THEN
[1]	101	!
	102	!-- Open file for time series output in NetCDF format
	103	dopts_time_count = dopts_time_count + 1
	104	CALL check_open( 109 )
	105	#if defined( __netcdf )
	106	!
	107	!-- Update the particle time series time axis
[291]	108	nc_stat = NF90_PUT_VAR( id_set_pts, id_var_time_pts, &
	109	(/ time_since_reference_point /), &
[1]	110	start = (/ dopts_time_count /), count = (/ 1 /) )
[263]	111	CALL handle_netcdf_error( 'data_output_ptseries', 391 )
[1]	112	#endif
	113
	114	ENDIF
	115
[825]	116	ALLOCATE( pts_value(0:number_of_particle_groups,dopts_num), &
	117	pts_value_l(0:number_of_particle_groups,dopts_num) )
	118
[1353]	119	pts_value_l = 0.0_wp
	120	pts_value_l(:,16) = 9999999.9_wp ! for calculation of minimum radius
[1]	121
	122	!
	123	!-- Calculate or collect the particle time series quantities for all particles
	124	!-- and seperately for each particle group (if there is more than one group)
[1359]	125	DO i = nxl, nxr
	126	DO j = nys, nyn
	127	DO k = nzb, nzt
	128	number_of_particles = prt_count(k,j,i)
	129	IF (number_of_particles <= 0) CYCLE
	130	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	131	DO n = 1, number_of_particles
[1]	132
[1359]	133	IF ( particles(n)%particle_mask ) THEN ! Restrict analysis to active particles
[1]	134
[1359]	135	pts_value_l(0,1) = pts_value_l(0,1) + 1.0_wp ! total # of particles
	136	pts_value_l(0,2) = pts_value_l(0,2) + &
	137	( particles(n)%x - particles(n)%origin_x ) ! mean x
	138	pts_value_l(0,3) = pts_value_l(0,3) + &
	139	( particles(n)%y - particles(n)%origin_y ) ! mean y
	140	pts_value_l(0,4) = pts_value_l(0,4) + &
	141	( particles(n)%z - particles(n)%origin_z ) ! mean z
	142	pts_value_l(0,5) = pts_value_l(0,5) + particles(n)%z ! mean z (absolute)
	143	pts_value_l(0,6) = pts_value_l(0,6) + particles(n)%speed_x ! mean u
	144	pts_value_l(0,7) = pts_value_l(0,7) + particles(n)%speed_y ! mean v
	145	pts_value_l(0,8) = pts_value_l(0,8) + particles(n)%speed_z ! mean w
	146	IF ( .NOT. curvature_solution_effects ) THEN
	147	pts_value_l(0,9) = pts_value_l(0,9) + particles(n)%rvar1 ! mean sgsu
	148	pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
	149	pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
	150	ENDIF
	151	IF ( particles(n)%speed_z > 0.0_wp ) THEN
	152	pts_value_l(0,12) = pts_value_l(0,12) + 1.0_wp ! # of upward moving prts
	153	pts_value_l(0,13) = pts_value_l(0,13) + &
	154	particles(n)%speed_z ! mean w upw.
	155	ELSE
	156	pts_value_l(0,14) = pts_value_l(0,14) + &
	157	particles(n)%speed_z ! mean w down
	158	ENDIF
	159	pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius ! mean rad
	160	pts_value_l(0,16) = MIN( pts_value_l(0,16), particles(n)%radius ) ! minrad
	161	pts_value_l(0,17) = MAX( pts_value_l(0,17), particles(n)%radius ) ! maxrad
	162	pts_value_l(0,18) = pts_value_l(0,18) + 1.0_wp
	163	pts_value_l(0,19) = pts_value_l(0,18) + 1.0_wp
[1]	164	!
[1359]	165	!-- Repeat the same for the respective particle group
	166	IF ( number_of_particle_groups > 1 ) THEN
	167	jg = particles(n)%group
[1]	168
[1359]	169	pts_value_l(jg,1) = pts_value_l(jg,1) + 1.0_wp
	170	pts_value_l(jg,2) = pts_value_l(jg,2) + &
	171	( particles(n)%x - particles(n)%origin_x )
	172	pts_value_l(jg,3) = pts_value_l(jg,3) + &
	173	( particles(n)%y - particles(n)%origin_y )
	174	pts_value_l(jg,4) = pts_value_l(jg,4) + &
	175	( particles(n)%z - particles(n)%origin_z )
	176	pts_value_l(jg,5) = pts_value_l(jg,5) + particles(n)%z
	177	pts_value_l(jg,6) = pts_value_l(jg,6) + particles(n)%speed_x
	178	pts_value_l(jg,7) = pts_value_l(jg,7) + particles(n)%speed_y
	179	pts_value_l(jg,8) = pts_value_l(jg,8) + particles(n)%speed_z
	180	IF ( .NOT. curvature_solution_effects ) THEN
	181	pts_value_l(jg,9) = pts_value_l(jg,9) + particles(n)%rvar1
	182	pts_value_l(jg,10) = pts_value_l(jg,10) + particles(n)%rvar2
	183	pts_value_l(jg,11) = pts_value_l(jg,11) + particles(n)%rvar3
	184	ENDIF
	185	IF ( particles(n)%speed_z > 0.0_wp ) THEN
	186	pts_value_l(jg,12) = pts_value_l(jg,12) + 1.0_wp
	187	pts_value_l(jg,13) = pts_value_l(jg,13) + particles(n)%speed_z
	188	ELSE
	189	pts_value_l(jg,14) = pts_value_l(jg,14) + particles(n)%speed_z
	190	ENDIF
	191	pts_value_l(jg,15) = pts_value_l(jg,15) + particles(n)%radius
	192	pts_value_l(jg,16) = MIN( pts_value(jg,16), particles(n)%radius )
	193	pts_value_l(jg,17) = MAX( pts_value(jg,17), particles(n)%radius )
	194	pts_value_l(jg,18) = pts_value_l(jg,18) + 1.0_wp
	195	pts_value_l(jg,19) = pts_value_l(jg,19) + 1.0_wp
	196	ENDIF
[1]	197
[1359]	198	ENDIF
[1]	199
[1359]	200	ENDDO
	201
	202	ENDDO
	203	ENDDO
[1]	204	ENDDO
	205
[1359]	206
[1]	207	#if defined( __parallel )
	208	!
	209	!-- Sum values of the subdomains
	210	inum = number_of_particle_groups + 1
	211
[622]	212	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[824]	213	CALL MPI_ALLREDUCE( pts_value_l(0,1), pts_value(0,1), 15*inum, MPI_REAL, &
[1]	214	MPI_SUM, comm2d, ierr )
[622]	215	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[824]	216	CALL MPI_ALLREDUCE( pts_value_l(0,16), pts_value(0,16), inum, MPI_REAL, &
[825]	217	MPI_MIN, comm2d, ierr )
	218	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	219	CALL MPI_ALLREDUCE( pts_value_l(0,17), pts_value(0,17), inum, MPI_REAL, &
[1]	220	MPI_MAX, comm2d, ierr )
[622]	221	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[825]	222	CALL MPI_ALLREDUCE( pts_value_l(0,18), pts_value(0,18), inum, MPI_REAL, &
	223	MPI_MAX, comm2d, ierr )
	224	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	225	CALL MPI_ALLREDUCE( pts_value_l(0,19), pts_value(0,19), inum, MPI_REAL, &
[1]	226	MPI_MIN, comm2d, ierr )
	227	#else
[825]	228	pts_value(:,1:19) = pts_value_l(:,1:19)
[1]	229	#endif
	230
	231	!
[825]	232	!-- Normalize the above calculated quantities (except min/max values) with the
	233	!-- total number of particles
[1]	234	IF ( number_of_particle_groups > 1 ) THEN
	235	inum = number_of_particle_groups
	236	ELSE
	237	inum = 0
	238	ENDIF
	239
	240	DO j = 0, inum
	241
[1353]	242	IF ( pts_value(j,1) > 0.0_wp ) THEN
[1]	243
[824]	244	pts_value(j,2:15) = pts_value(j,2:15) / pts_value(j,1)
[1353]	245	IF ( pts_value(j,12) > 0.0_wp .AND. pts_value(j,12) < 1.0_wp ) THEN
[1]	246	pts_value(j,13) = pts_value(j,13) / pts_value(j,12)
[1353]	247	pts_value(j,14) = pts_value(j,14) / ( 1.0_wp - pts_value(j,12) )
	248	ELSEIF ( pts_value(j,12) == 0.0_wp ) THEN
	249	pts_value(j,13) = -1.0_wp
[1]	250	ELSE
[1353]	251	pts_value(j,14) = -1.0_wp
[1]	252	ENDIF
	253
	254	ENDIF
	255
	256	ENDDO
	257
	258	!
	259	!-- Calculate higher order moments of particle time series quantities,
	260	!-- seperately for each particle group (if there is more than one group)
[1359]	261	DO i = nxl, nxr
	262	DO j = nys, nyn
	263	DO k = nzb, nzt
	264	number_of_particles = prt_count(k,j,i)
	265	IF (number_of_particles <= 0) CYCLE
	266	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	267	DO n = 1, number_of_particles
[1]	268
[1359]	269	pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%x - &
	270	particles(n)%origin_x - pts_value(0,2) )*2 ! x2
	271	pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%y - &
	272	particles(n)%origin_y - pts_value(0,3) )*2 ! y2
	273	pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%z - &
	274	particles(n)%origin_z - pts_value(0,4) )*2 ! z2
	275	pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_x - &
	276	pts_value(0,6) )*2 ! u2
	277	pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%speed_y - &
	278	pts_value(0,7) )*2 ! v2
	279	pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - &
	280	pts_value(0,8) )*2 ! w2
	281	IF ( .NOT. curvature_solution_effects ) THEN
	282	pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
	283	pts_value(0,9) )**2 ! u"2
	284	pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
	285	pts_value(0,10) )**2 ! v"2
	286	pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
	287	pts_value(0,11) )**2 ! w"2
	288	ENDIF
[1]	289	!
[1359]	290	!-- Repeat the same for the respective particle group
	291	IF ( number_of_particle_groups > 1 ) THEN
	292	jg = particles(n)%group
[1]	293
[1359]	294	pts_value_l(jg,20) = pts_value_l(jg,20) + ( particles(n)%x - &
	295	particles(n)%origin_x - pts_value(jg,2) )**2
	296	pts_value_l(jg,21) = pts_value_l(jg,21) + ( particles(n)%y - &
	297	particles(n)%origin_y - pts_value(jg,3) )**2
	298	pts_value_l(jg,22) = pts_value_l(jg,22) + ( particles(n)%z - &
	299	particles(n)%origin_z - pts_value(jg,4) )**2
	300	pts_value_l(jg,23) = pts_value_l(jg,23) + ( particles(n)%speed_x - &
	301	pts_value(jg,6) )**2
	302	pts_value_l(jg,24) = pts_value_l(jg,24) + ( particles(n)%speed_y - &
	303	pts_value(jg,7) )**2
	304	pts_value_l(jg,25) = pts_value_l(jg,25) + ( particles(n)%speed_z - &
	305	pts_value(jg,8) )**2
	306	IF ( .NOT. curvature_solution_effects ) THEN
	307	pts_value_l(jg,26) = pts_value_l(jg,26) + ( particles(n)%rvar1 - &
	308	pts_value(jg,9) )**2
	309	pts_value_l(jg,27) = pts_value_l(jg,27) + ( particles(n)%rvar2 - &
	310	pts_value(jg,10) )**2
	311	pts_value_l(jg,28) = pts_value_l(jg,28) + ( particles(n)%rvar3 - &
	312	pts_value(jg,11) )**2
	313	ENDIF
	314	ENDIF
[1]	315
[1359]	316	ENDDO
	317	ENDDO
	318	ENDDO
[1]	319	ENDDO
	320
[825]	321	pts_value_l(0,29) = ( number_of_particles - pts_value(0,1) / numprocs )**2
[1]	322	! variance of particle numbers
	323	IF ( number_of_particle_groups > 1 ) THEN
	324	DO j = 1, number_of_particle_groups
[825]	325	pts_value_l(j,29) = ( pts_value_l(j,1) - &
[1]	326	pts_value(j,1) / numprocs )**2
	327	ENDDO
	328	ENDIF
	329
	330	#if defined( __parallel )
	331	!
	332	!-- Sum values of the subdomains
	333	inum = number_of_particle_groups + 1
	334
[622]	335	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[825]	336	CALL MPI_ALLREDUCE( pts_value_l(0,20), pts_value(0,20), inum*10, MPI_REAL, &
[1]	337	MPI_SUM, comm2d, ierr )
	338	#else
[825]	339	pts_value(:,20:29) = pts_value_l(:,20:29)
[1]	340	#endif
	341
	342	!
	343	!-- Normalize the above calculated quantities with the total number of
	344	!-- particles
	345	IF ( number_of_particle_groups > 1 ) THEN
	346	inum = number_of_particle_groups
	347	ELSE
	348	inum = 0
	349	ENDIF
	350
	351	DO j = 0, inum
	352
[1353]	353	IF ( pts_value(j,1) > 0.0_wp ) THEN
[825]	354	pts_value(j,20:28) = pts_value(j,20:28) / pts_value(j,1)
[1]	355	ENDIF
[825]	356	pts_value(j,29) = pts_value(j,29) / numprocs
[1]	357
	358	ENDDO
	359
	360	#if defined( __netcdf )
	361	!
	362	!-- Output particle time series quantities in NetCDF format
[1327]	363	IF ( myid == 0 ) THEN
[1]	364	DO j = 0, inum
	365	DO i = 1, dopts_num
	366	nc_stat = NF90_PUT_VAR( id_set_pts, id_var_dopts(i,j), &
	367	(/ pts_value(j,i) /), &
	368	start = (/ dopts_time_count /), &
	369	count = (/ 1 /) )
[263]	370	CALL handle_netcdf_error( 'data_output_ptseries', 392 )
[1]	371	ENDDO
	372	ENDDO
	373	ENDIF
	374	#endif
	375
[825]	376	DEALLOCATE( pts_value, pts_value_l )
	377
[1318]	378	CALL cpu_log( log_point(36), 'data_output_ptseries', 'stop' )
[1]	379
	380	END SUBROUTINE data_output_ptseries

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |