Home

Context Navigation

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

Last change on this file since 1320 was 1320, checked in by raasch, 11 years ago
ONLY-attribute added to USE-statements, kind-parameters added to all INTEGER and REAL declaration statements, kinds are defined in new module kinds, old module precision_kind is removed, revision history before 2012 removed, comment fields (!:) to be used for variable explanations added to all variable declaration statements
Property svn:keywords set to `Id`
File size: 14.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |