source: palm/trunk/UTIL/read_prt_data.f90 @ 2125

Last change on this file since 2125 was 2125, checked in by hoffmann, 8 years ago

last commit documented

File size: 6.9 KB
Line 
1 PROGRAM read_prt_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-2014  Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: read_prt_data.f90 2123 2017-01-18 12:49:59Z hoffmann $
26!
27! 2123 2017-01-18 12:49:59Z hoffmann
28! Particle variables renamed: tailpoints, tail_id, and dvrp_psize to id1, id2,
29! and user
30!
31! 1771 2016-02-29 10:57:56Z hoffmann
32! Initial revision.
33!
34! Description:
35! ------------
36! This routine reads the particle data generated by PALM, and enables user
37! analysis. Compile and execute this program in the folder where your particle
38! data (_000000, _000001, ...) is stored.
39!-------------------------------------------------------------------------------!
40
41    IMPLICIT NONE
42
43    CHARACTER(LEN=7)    ::  i_proc_char
44    CHARACTER (LEN=80)  ::  rtext
45    CHARACTER (LEN=110) ::  run_description_header
46
47    REAL(KIND=8) ::  simulated_time
48    INTEGER      ::  ip, i_proc=0, i_proc_end, jp, kp,                         &
49                     max_number_of_particle_groups, number_of_particles,       &
50                     number_of_particle_groups, n, nxl,                        &
51                     nxr, nys, nyn, nzb, nzt, nbgp, status
52
53    INTEGER, DIMENSION(:,:,:), ALLOCATABLE ::  prt_count
54
55    TYPE particle_type
56       SEQUENCE
57       REAL(KIND=8) ::  radius, age, age_m, dt_sum, user, e_m,                 &
58                        origin_x, origin_y, origin_z, rvar1, rvar2, rvar3,     &
59                        speed_x, speed_y, speed_z, weight_factor, x, y, z
60      INTEGER       ::  class, group, id1, id2
61      LOGICAL       ::  particle_mask
62      INTEGER       ::  block_nr
63    END TYPE particle_type
64
65    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  particles_temp
66    TYPE(particle_type), DIMENSION(:), POINTER     ::  particles
67
68    TYPE  grid_particle_def
69       TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  particles
70    END TYPE grid_particle_def
71
72    TYPE(grid_particle_def), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  grid_particles
73
74    TYPE particle_groups_type
75        SEQUENCE
76        REAL(KIND=8) ::  density_ratio, radius, exp_arg, exp_term
77    END TYPE particle_groups_type
78
79    TYPE(particle_groups_type), DIMENSION(:), ALLOCATABLE ::  particle_groups
80
81    LOGICAL ::  found
82!
83!-- Check if file from PE0 exists and terminate program if it doesn't.
84    WRITE (i_proc_char,'(''_'',I6.6)')  i_proc
85    INQUIRE ( FILE=i_proc_char, EXIST=found )
86!
87!-- Estimate the number of files (equal to the number of PEs which
88!-- have been used in PALM)
89    DO  WHILE ( found )
90       OPEN ( 85, FILE=i_proc_char, FORM='UNFORMATTED' )
91       i_proc = i_proc + 1
92       WRITE (i_proc_char,'(''_'',I6.6)')  i_proc
93       INQUIRE ( FILE=i_proc_char, EXIST=found )
94       CLOSE( 85 )
95    ENDDO
96!
97!-- Info-output
98    PRINT*, ''
99    PRINT*, '*** read_prt ***'
100    IF ( i_proc /= 0 )  THEN
101       PRINT*, '***', i_proc, ' file(s) found'
102    ELSE
103       PRINT*, '+++ file _000000 not found'
104       PRINT*, '    program terminated'
105       STOP
106    ENDIF
107!
108!-- Set number of files and opens them sequentially
109    i_proc_end = i_proc-1
110
111    DO  i_proc = 0, i_proc_end
112!
113!--    Open particle data file for each process
114       WRITE (i_proc_char,'(''_'',I6.6)')  i_proc
115       OPEN ( 85, FILE=i_proc_char, FORM='UNFORMATTED' ) 
116!
117!--    Read general description of file and inform user
118       READ ( 85 )  run_description_header
119       READ ( 85 )  rtext
120       IF ( i_proc == 0 )  THEN
121          PRINT*, ' '
122          PRINT*, '*** data description header:'
123          PRINT*, '*** ', run_description_header
124          PRINT*, '*** ', rtext
125          PRINT*, ' '
126       ENDIF
127       PRINT*, '*** data of file', i_proc, 'is analysed'
128       READ ( 85 )  number_of_particle_groups, max_number_of_particle_groups
129
130       IF ( .NOT. ALLOCATED(particle_groups) )  THEN
131          ALLOCATE( particle_groups(max_number_of_particle_groups) )
132       ENDIF
133   
134       READ ( 85 )  particle_groups
135       READ ( 85 )  nxl, nxr, nys, nyn, nzb, nzt, nbgp
136
137       IF ( ALLOCATED(prt_count)  .OR.  ALLOCATED(grid_particles) )  THEN
138          DEALLOCATE( prt_count, grid_particles )
139       ENDIF
140
141       ALLOCATE( prt_count(nzb:nzt+1,nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp),     &
142                 grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )
143
144!
145!--    Start individual time loop
146       DO
147          READ( 85, IOSTAT=status)  simulated_time
148          IF ( status < 0 )  THEN
149             PRINT*, '*** end of file reached'
150             EXIT
151          ENDIF
152          PRINT*, '*** time of analyzed data set:', simulated_time
153
154          READ ( 85 )  prt_count
155
156          DO  ip = nxl, nxr
157             DO  jp = nys, nyn
158                DO  kp = nzb+1, nzt
159                   number_of_particles = prt_count(kp,jp,ip)
160                   IF ( number_of_particles <= 0 )  CYCLE
161                   ALLOCATE( grid_particles(kp,jp,ip)%particles(1:number_of_particles) )
162                   ALLOCATE( particles_temp(1:number_of_particles) )
163                   READ ( 85 )  particles_temp
164                   grid_particles(kp,jp,ip)%particles = particles_temp
165                   DEALLOCATE( particles_temp )
166                ENDDO
167             ENDDO
168          ENDDO
169!
170!--       This part can be used for user analysis
171          DO  ip = nxl, nxr
172             DO  jp = nys, nyn
173                DO  kp = nzb+1, nzt
174                   number_of_particles = prt_count(kp,jp,ip)
175                   IF ( number_of_particles <= 0 )  CYCLE
176                   particles => grid_particles(kp,jp,ip)%particles
177!
178!--                Add your analysis here
179!                   DO  n = 1, number_of_particles
180!
181!                   ENDDO
182
183                ENDDO
184             ENDDO
185          ENDDO
186!
187!--       Deallocate grid_particles%particles in case of more than one timestep
188          DO  ip = nxl, nxr
189             DO  jp = nys, nyn
190                DO  kp = nzb+1, nzt
191                   number_of_particles = prt_count(kp,jp,ip)
192                   IF ( number_of_particles <= 0 )  CYCLE
193                   DEALLOCATE( grid_particles(kp,jp,ip)%particles )
194                ENDDO
195             ENDDO
196          ENDDO     
197
198       ENDDO ! end of time loop
199
200       CLOSE( 85 )
201
202    ENDDO ! end of file loop
203
204 END PROGRAM read_prt_data
Note: See TracBrowser for help on using the repository browser.