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

Last change on this file since 4010 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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