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

Last change on this file since 2065 was 1772, checked in by hoffmann, 9 years ago

last commit documented

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