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

Last change on this file since 3485 was 2978, checked in by schwenkel, 6 years ago

Bugfix in particle analyze tool

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