Home

Context Navigation

data_output_ptseries.f90 @ 1327

Last change on this file since 1327 was 1327, checked in by raasch, 11 years ago

Changed:

-s real64 removed (.mrun.config.hlrnIII)
-r8 removed (.mrun.config.imuk)
deleted: .mrun.config.imuk_ice2_netcdf4 .mrun.config.imuk_hlrn

REAL constants defined as wp-kind in modules

"baroclinicity" renamed "baroclinity", "ocean version" replaced by
"ocean mode"

code parts concerning old output formats "iso2d" and "avs" removed.
netCDF is the only remaining output format.

Errors:

bugfix: duplicate error message 56 removed

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

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

Context Navigation

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

Download in other formats: