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

Last change on this file since 635 was 1, checked in by raasch, 18 years ago

Initial repository layout and content

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