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

Last change on this file since 1046 was 1046, checked in by maronga, 12 years ago

put scripts and utilities under GPL

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