source: palm/trunk/UTIL/analyze_particle_data.f90

Last change on this file was 4843, checked in by raasch, 23 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

  • 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-2021  Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: analyze_particle_data.f90 4843 2021-01-15 15:22:11Z banzhafs $
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.