source: palm/trunk/UTIL/analyze_particle_data.f90 @ 1933

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

update of GPL copyright

  • Property svn:keywords set to Id
File size: 10.4 KB
Line 
1 PROGRAM analyze_particle_data
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2014  Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22! Former revisions:
23! -----------------
24! $Id: analyze_particle_data.f90 1310 2014-03-14 08:01:56Z hellstea $
25!
26! 1046 2012-11-09 14:38:45Z maronga
27! code put under GPL (PALM 3.9)
28!
29! Description:
30! ------------
31! This routine reads the particle data files generated by PALM
32! and does some statistical analysis on these data.
33!-------------------------------------------------------------------------------!
34
35    IMPLICIT NONE
36
37!
38!-- Variable definitions
39    CHARACTER (LEN=5)  ::  id_char
40    CHARACTER (LEN=80) ::  run_description_header
41
42    INTEGER, PARAMETER ::  spk = SELECTED_REAL_KIND( 6 )
43
44    INTEGER ::  class, danz = 0, i, id, maximum_number_of_particles, n,   &
45                number_of_intervalls, number_of_particles,                &
46                number_of_vertical_levels, timelevel = 1, vertical_level, &
47                total_number_of_particles
48
49    INTEGER, DIMENSION(:), ALLOCATABLE ::  class_table
50   
51    LOGICAL ::  found
52
53    REAL    ::  class_width, hdistance, km, km_x, km_y, particle_age, sigma, &
54                sigma_local, sigma_local_x, sigma_local_y, sigma_x, sigma_y, &
55                simulated_time, vertical_resolution
56    REAL, DIMENSION(:,:), ALLOCATABLE ::  diffusivities
57
58    TYPE particle_type
59       SEQUENCE
60       INTEGER ::  color, tailpoints
61       REAL    ::  age, origin_x, origin_y, origin_z, size, speed_x, speed_y, &
62                   speed_z, x, y, z
63    END TYPE particle_type
64
65    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  particles
66
67
68!
69!-- Check, if file from PE0 exists and terminate program if it doesn't.
70    WRITE (id_char,'(''_'',I4.4)')  danz
71    INQUIRE ( FILE=id_char, EXIST=found )
72!
73!-- Find out the number of files (equal to the number of PEs which
74!-- have been used in PALM) and open them
75    DO  WHILE ( found )
76
77       OPEN ( danz+110, FILE=id_char, FORM='UNFORMATTED' )
78       danz = danz + 1
79       WRITE (id_char,'(''_'',I4.4)')  danz
80       INQUIRE ( FILE=id_char, EXIST=found )
81
82    ENDDO
83!
84!-- Info-output
85    PRINT*, ''
86    PRINT*, '*** analyze_particle_data ***'
87    IF ( danz /= 0 )  THEN
88       PRINT*, '*** ', danz, ' file(s) found'
89    ELSE
90       PRINT*, '+++ file _0000 not found'
91       PRINT*, '    program terminated'
92       STOP
93    ENDIF
94
95!
96!-- Loop over all timelevels of output
97    DO
98       total_number_of_particles = 0
99       sigma   = 0.0
100       sigma_x = 0.0
101       sigma_y = 0.0
102!
103!--    Loop over all files (reading data of the subdomains)
104       DO  id = 0, danz-1
105!
106!--       Read file header
107          IF ( timelevel == 1 )  THEN
108             READ ( id+110 )  run_description_header
109!
110!--          Print header information
111             IF ( id == 0 )  THEN
112                PRINT*, '*** run: ', run_description_header
113                PRINT*, ' '
114                PRINT*, '--> enter class width in m:'
115                READ*, class_width
116                PRINT*, '--> enter number of class intervalls:'
117                READ*, number_of_intervalls
118                PRINT*, '--> enter vertical resolution in m:'
119                READ*, vertical_resolution
120                PRINT*, '--> enter number of vertical levels:'
121                READ*, number_of_vertical_levels
122!
123!--             Allocate table space
124                ALLOCATE( class_table( 0:number_of_intervalls ) )
125                class_table = 0
126                ALLOCATE( diffusivities(0:number_of_vertical_levels,5) )
127                diffusivities = 0.0
128             ENDIF
129          ENDIF
130!
131!--       Read time information and indices
132          READ ( id+110, END=10 )  simulated_time, maximum_number_of_particles,&
133                                   number_of_particles
134!
135!--       Print timelevel and number of particles
136          IF ( id == 0 )  THEN
137             PRINT*, ' '
138             PRINT*, '*** time: ', simulated_time
139          ENDIF
140          PRINT*, 'PE', id, ': ', number_of_particles, ' particles'
141!
142!--       Allocate array and read particle data
143          ALLOCATE( particles(maximum_number_of_particles) )
144          READ ( id+110 )  particles
145!
146!--       Analyze the particle data
147          DO  n = 1, number_of_particles
148!
149!--          Calculate horizontal distance from of particle from its origin
150             hdistance = SQRT( ( particles(n)%x - particles(n)%origin_x )**2 + &
151                               ( particles(n)%y - particles(n)%origin_y )**2 )
152             class = hdistance / class_width
153             sigma_local   = hdistance**2
154             sigma_local_x = ( particles(n)%x - particles(n)%origin_x )**2
155             sigma_local_y = ( particles(n)%y - particles(n)%origin_y )**2
156
157             vertical_level = particles(n)%origin_z / vertical_resolution
158             IF ( vertical_level > number_of_vertical_levels )  THEN
159                vertical_level = number_of_vertical_levels
160             ENDIF
161
162             IF ( class > number_of_intervalls )  THEN
163                class = number_of_intervalls
164!                PRINT*, 'x =',particles(n)%x,' y =',particles(n)%y
165!                PRINT*, 'xo=',particles(n)%origin_x,' yo=',particles(n)%origin_y
166             ENDIF
167
168             class_table(class) = class_table(class) + 1
169
170             diffusivities(vertical_level,1) = diffusivities(vertical_level,1) +&
171                                               sigma_local
172             diffusivities(vertical_level,2) = diffusivities(vertical_level,2) +&
173                                               sigma_local_x
174             diffusivities(vertical_level,3) = diffusivities(vertical_level,3) +&
175                                               sigma_local_y
176             diffusivities(vertical_level,4) = diffusivities(vertical_level,4) +&
177                                               1.0
178
179             vertical_level = particles(n)%z / vertical_resolution
180             IF ( vertical_level > number_of_vertical_levels )  THEN
181                vertical_level = number_of_vertical_levels
182             ENDIF
183             diffusivities(vertical_level,5) = diffusivities(vertical_level,5) +&
184                                               1.0
185
186!
187!--          Summation for variances
188             sigma = sigma + sigma_local
189             sigma_x = sigma_x + sigma_local_x
190             sigma_y = sigma_y + sigma_local_y
191             total_number_of_particles = total_number_of_particles + 1
192
193          ENDDO
194!
195!--       Store the particle age (it is provided that all particles have the
196!--       same age)
197          particle_age = particles(1)%age
198
199!
200!--       Deallocate particle array before data from next file are read
201          DEALLOCATE( particles )
202
203       ENDDO  ! next file
204!
205!--    Print statistics
206       PRINT*, ' '
207       PRINT*, '*** statistics for t = ', simulated_time
208       DO  n = 0, number_of_intervalls-1
209          WRITE ( *, 1 ) n*class_width, (n+1)*class_width, class_table(n)
210   1      FORMAT (F6.1,' - ',F6.1, ' m   n = ',I7)
211       ENDDO
212       WRITE ( *, 2 )  (number_of_intervalls+1)*class_width, &
213                       class_table(number_of_intervalls)
214   2   FORMAT (6X,' > ',F6.1,' m   n = ',I7)
215
216       sigma   = SQRT( sigma   / REAL( total_number_of_particles ) )
217       km      = sigma**2 / ( 2.0 * particle_age )
218       sigma_x = SQRT( sigma_x / REAL( total_number_of_particles ) )
219       km_x    = sigma_x**2 / ( 2.0 * particle_age )
220       sigma_y = SQRT( sigma_y / REAL( total_number_of_particles ) )
221       km_y    = sigma_y**2 / ( 2.0 * particle_age )
222       PRINT*, ' '
223       WRITE ( *, 3 )  sigma, km, sigma_x, km_x, sigma_y, km_y
224   3   FORMAT ('sigma   = ',F6.1,' m   Km   = ',F5.1,' m**2/s'/ &
225               'sigma_x = ',F6.1,' m   Km_x = ',F5.1,' m**2/s'/ &
226               'sigma_y = ',F6.1,' m   Km_y = ',F5.1,' m**2/s')
227
228       PRINT*, ' '
229       PRINT*, 'Height dependence of diffusivities:'
230       DO  i = 0, number_of_vertical_levels-1
231          IF ( diffusivities(i,4) == 0.0 )  diffusivities(i,4) = 1.0E-20
232          WRITE ( *, 4 )  i*vertical_resolution, (i+1.0)*vertical_resolution,&
233                          ( diffusivities(i,1) / diffusivities(i,4) ) / &
234                                     ( 2.0 * particle_age ),            &
235                          ( diffusivities(i,2) / diffusivities(i,4) ) / &
236                                     ( 2.0 * particle_age ),            &
237                          ( diffusivities(i,3) / diffusivities(i,4) ) / &
238                                     ( 2.0 * particle_age ),            &
239                          diffusivities(i,4), diffusivities(i,5)
240   4      FORMAT (F6.1,'-',F6.1,' m  Km=',F5.1,' Km_x=',F5.1, &
241                  ' Km_y=',F5.1,' n_o=',F7.0,' n=',F7.0)
242       ENDDO
243       IF ( diffusivities(number_of_vertical_levels,4) == 0.0 )  THEN
244          diffusivities(number_of_vertical_levels,4) = 1.0E-20
245       ENDIF
246          i = number_of_vertical_levels
247          WRITE ( *, 5 )  i*vertical_resolution,                        &
248                          ( diffusivities(i,1) / diffusivities(i,4) ) / &
249                                     ( 2.0 * particle_age ),            &
250                          ( diffusivities(i,2) / diffusivities(i,4) ) / &
251                                     ( 2.0 * particle_age ),            &
252                          ( diffusivities(i,3) / diffusivities(i,4) ) / &
253                                     ( 2.0 * particle_age ),            &
254                          diffusivities(i,4), diffusivities(i,5)
255   5   FORMAT (F6.1,'-...... m  Km=',F5.1,' Km_x=',F5.1, &
256                  ' Km_y=',F5.1,' n_o=',F7.0,' n=',F7.0)
257
258!
259!--    Initialize class table for next timelevel
260       class_table   =   0
261       diffusivities = 0.0
262       timelevel = timelevel + 1
263       
264    ENDDO  ! next timelevel
265
26610  PRINT*, '*** EOF reached on file PARTICLE_DATA/_0000'
267
268 END PROGRAM analyze_particle_data
Note: See TracBrowser for help on using the repository browser.