Home

Context Navigation

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

Last change on this file since 947 was 826, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 12.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |