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

Last change on this file since 2698 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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