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

Last change on this file since 4417 was 4262, checked in by schwenkel, 4 years ago

Added BIND attribute and C datatypes

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