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

Last change on this file since 368 was 291, checked in by raasch, 16 years ago

changes for coupling with independent precursor runs; z_i calculation with Sullivan criterion

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