source: palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90 @ 3496

Last change on this file since 3496 was 3496, checked in by suehring, 5 years ago

Use subroutine call for intrinsic routine fseek instead of function call. gfortran has some problems with the function interface.

  • Property svn:keywords set to Id
File size: 17.4 KB
Line 
1!> @file surface_output_merge.f90
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: surface_output_to_vtk.f90 3496 2018-11-06 15:59:50Z suehring $
27! Use subroutine call for fseek instead of function call. gfortran has some
28! problems with this.
29!
30! 3494 2018-11-06 14:51:27Z suehring
31! Initial version
32!
33! Authors:
34! --------
35! @author Matthias Suehring and Klaus Ketelsen
36!
37!------------------------------------------------------------------------------!
38! Description:
39! ------------
40!> This routine combines surface output from PALM-subdomains into one file.
41!> Output from every processor element is opened and read and output
42!> from all processor elements are written into one file for each timestep
43!> according to Paraview's VTK standard.
44!> Output is distinguished between instantaneous and time-averaged data.
45!------------------------------------------------------------------------------!
46 PROGRAM surface_output_to_vtk
47
48    IMPLICIT NONE
49   
50    CHARACTER(LEN=4)   ::  char_time            !< string indicating simulated time
51    CHARACTER(LEN=4)   ::  file_suffix = '.bin' !< string which contain the suffix indicating surface data
52   
53    CHARACTER(LEN=10)  ::  char_dum             !< dummy string
54   
55    CHARACTER(LEN=30)  ::  myid_char            !< combined string indicating binary file
56   
57    CHARACTER(LEN=100) ::  path                 !< path to the binary data
58    CHARACTER(LEN=100) ::  run                  !< name of the run
59    CHARACTER(LEN=100) ::  variable_name        !< name of the processed output variable   
60   
61    INTEGER(4)   ::  ftell                      !< intrinsic function, get current position in file
62!     INTEGER(4)   ::  fseek                      !< intrinsic function, go to given position in file
63    INTEGER(4)   ::  ndum                       !< return parameter of intrinsic function fseek
64   
65    INTEGER, PARAMETER ::  iwp = 4              !< integer precision
66    INTEGER, PARAMETER ::  wp  = 8              !< float precision
67   
68    INTEGER(iwp) ::  cycle_number               !< cycle number
69    INTEGER(iwp) ::  f                          !< running index over all binary files
70    INTEGER(iwp) ::  file_id_in = 18            !< file unit for input binaray file   
71    INTEGER(iwp) ::  file_id_out = 20           !< file unit for output VTK file       
72    INTEGER(iwp) ::  file_id_out_header = 19    !< file unit for temporary header file
73    INTEGER(iwp) ::  length                     !< length of string on file
74    INTEGER(iwp) ::  n                          !< running index over all points and polygons
75    INTEGER(iwp) ::  npoints_total              !< total number of points
76    INTEGER(iwp) ::  ns_total                   !< total number of polygons
77    INTEGER(iwp) ::  num_pe                     !< number of processors of the run
78   
79    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  filepos !< current fileposition in binary file
80    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  npoints !< number of points/vertices in a binaray file
81    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns      !< number of surface elements in a binaray file
82   
83    LOGICAL ::  convert_average_data = .FALSE.          !< namelist parameter to decide whether average or instantaneous data should be converted
84    LOGICAL, DIMENSION(:), ALLOCATABLE      ::  eof     !< flag indicating that end of binary file is reached
85   
86    REAL(wp)                              ::  simulated_time !< output time   
87   
88    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var            !< actual surface data
89   
90    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  points         !< point / vertex data
91    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  polygons       !< polygon data
92   
93    logical :: flag
94
95!
96!-- Read namelist.
97    CALL surface_output_parin
98!
99!-- Allocate array which contains the file position in each output file,
100!-- in order to skip some reading.
101    ALLOCATE( eof( 0:num_pe-1 )     )
102    ALLOCATE( filepos( 0:num_pe-1 ) )
103    ALLOCATE( npoints( 0:num_pe-1 ) )
104    ALLOCATE( ns( 0:num_pe-1 )      )
105!
106!-- Initialize file position.
107    filepos = 0
108!
109!-- Open a temporary file which contains all necessary header information for the
110!-- VTK format.
111    OPEN ( file_id_out_header, FILE = 'HEADER', STATUS = 'REPLACE',            &
112           FORM = 'FORMATTED' ) 
113!
114!-- READ grid setup, i.e. the number and position of vertices and surface elements
115!-- and merge all this information into one file. Further, create all the required
116!-- header information of the VTK file.
117!-- Note, PARAVIEW expects one VTK file for each variable and each timestep.
118!-- Hence, header information needs to be duplicated multiple times, which will be
119!-- be done later in a bash script.
120!-- Moreover, Paraview expects consecutive vertex and polygon data, which are
121!-- all distributed over the binaray files. Hence, first read vertex data from
122!-- binary file, write this to the HEADER file, close the binary file and read
123!-- data from the next binary file and so on. This requires several openenings
124!-- and closings of the binary files and temporarily storage of the
125!-- file-positions.
126    DO  f = 0, num_pe - 1
127!
128!--    Create filename of the treated binary file.
129       CALL surface_output_create_file_string
130!
131!--    Open file with surface output for processor f.
132       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
133              TRIM( myid_char ), FORM = 'UNFORMATTED' )
134!
135!--    Read number of vertices / points and surface elements
136       READ ( file_id_in )  npoints(f)
137       READ ( file_id_in )  npoints_total
138       READ ( file_id_in )  ns(f)
139       READ ( file_id_in )  ns_total
140
141!
142!--    Allocate arrays where all the surfaces and vertices will be stored.
143       ALLOCATE( points(3,1:npoints(f))   )
144       
145!
146!--    Read polygon data and store them in a temporary file.
147       READ ( file_id_in )  points
148!
149!--    Obtain current file position. Will be stored for next file opening.
150       filepos(f) = ftell( file_id_in )
151!        CALL FTELL( file_id_in, filepos(f) )
152!
153!--    Write header information. Only one time required.
154       IF ( f == 0 )  THEN
155          WRITE( file_id_out_header,'(A)' ) "# vtk DataFile Version 3.0"
156          WRITE( file_id_out_header,'(A,F8.2,A)' ) "legacy vtk File generated by PALM,  simulation time = xxx sec"
157          WRITE( file_id_out_header,'(A)') "ASCII"
158
159          WRITE( file_id_out_header,'(A)') "DATASET POLYDATA"
160          WRITE( file_id_out_header,'(A,I5,A)') "POINTS ", npoints_total, " float"
161       ENDIF     
162!
163!--    Write the vertex data into header file.
164       DO  n = 1, npoints(f)
165          WRITE( file_id_out_header, '(8F10.1)' )  points(1:3,n)
166       ENDDO
167!
168!--    Deallocate vertex data and close binary file.
169       DEALLOCATE( points )
170       
171       CLOSE ( file_id_in )
172    ENDDO
173!
174!-- Now, treat polygon data.
175    DO  f = 0, num_pe - 1
176!
177!--    Create filename of the treated binary file .   
178       CALL surface_output_create_file_string
179!
180!--    Open file with surface output for processor f.
181       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
182              TRIM( myid_char ), FORM = 'UNFORMATTED' )
183!
184!--    Move to last postion.
185!        ndum = fseek( file_id_in, filepos(f), 0 )
186       CALL FSEEK( file_id_in, filepos(f), 0, ndum )
187!
188!--    Allocate array for polygon data
189       ALLOCATE( polygons(5,1:ns(f)) )
190!
191!--    Read polygon data and store them in a temporary file.
192       READ ( file_id_in )  polygons
193!
194!--    Obtain current file position after reading the local polygon data.
195!--    Will be used for next file opening.
196       filepos(f) = ftell( file_id_in )
197!        CALL FTELL( file_id_in, filepos(f) )
198!
199!--    Write further header information. Only one time required.
200       IF ( f == 0 )                                                           &
201          WRITE ( file_id_out_header, '(A,8I10)') "POLYGONS ",                 &
202                                                  ns_total, 5 * ns_total
203!
204!--    Write the polygons into the header file. Note, polygon data is of type
205!--    integer, as it just connects the point-indices which describe the given
206!--    surface element.
207       DO n = 1, ns(f)
208          WRITE ( file_id_out_header, '(8I10)' )  INT( polygons(1:5,n) )
209       ENDDO
210!
211!--    Deallocate array for polygon data and close the file.
212       DEALLOCATE( polygons )
213       
214       CLOSE ( file_id_in )
215       
216    ENDDO
217   
218    f = 0
219    CALL surface_output_create_file_string
220!
221!-- Write further header information. Only once required.
222    WRITE ( file_id_out_header, '(A,I10)') "CELL_DATA ", ns_total
223    WRITE ( file_id_out_header, '(A,I10)') "SCALARS cell_scalars float 1 "
224    WRITE ( file_id_out_header, '(A,I10)') "LOOKUP_TABLE default "
225!
226!-- Header creation has now been completed. File can be closed.   
227    CLOSE( file_id_out_header )
228!
229!-- Now, read all the actual data from the surface output. Please note, Paraview
230!-- VTK format expects 1 file per variable and time step.
231!-- The output file is only created once and includes the variable and the
232!-- simulated time.
233!-- In the binaray files, several variables and timesteps are stored, data for a
234!-- given variable, however, is distributed over all binary files. Hence, read
235!-- variable data for a given timestep from the binary file, write this data into
236!-- the target output file, remember file position in binary file and close it,
237!-- open nex binary file and read the variable, and so on, until all variables
238!-- for all timesteps are processed.
239    eof = .FALSE.
240    DO
241       DO  f = 0, num_pe - 1
242!
243!--       Clean up strings-
244          char_time = ''
245          variable_name = ''
246!
247!--       Create filename of the treated binary file.           
248          CALL surface_output_create_file_string
249!
250!--       Open binary file with surface output for processor f.
251          OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //         &
252                 TRIM( myid_char ), FORM = 'UNFORMATTED' )
253!
254!--       Move to last postion.
255!           ndum = fseek( file_id_in, filepos(f), 0 )
256          CALL FSEEK( file_id_in, filepos(f), 0, ndum )
257!
258!--       Read string length and string indicating the output time step.
259          READ ( file_id_in ) length
260          READ ( file_id_in ) char_time(1:length)
261!
262!--       If string for the output time indicates the end-of-file, set the eof
263!--       flag and skip the read of the loop.
264          IF ( char_time(1:length) == 'END' )  THEN
265             eof(f) = .TRUE. 
266             CLOSE ( file_id_in )
267             CYCLE
268          ENDIF
269!
270!--       Read output time, and variable name.
271          READ ( file_id_in ) simulated_time
272          READ ( file_id_in ) length
273          READ ( file_id_in ) variable_name(1:length)
274!
275!--       For first loop index, open the target output file. First create the
276!--       filename string. Further, copy HEADER file with the given filename
277!--       string. The header information must be given in each VTK file!
278          IF ( f == 0 )  THEN
279             IF ( simulated_time < 10.0_wp )  THEN
280                WRITE( char_dum, '(I1)' )  INT( simulated_time ) 
281             ELSEIF ( simulated_time < 100.0_wp )  THEN
282                WRITE( char_dum, '(I2)' )  INT( simulated_time )
283             ELSEIF ( simulated_time < 1000.0_wp )  THEN
284                WRITE( char_dum, '(I3)' )  INT( simulated_time )
285             ELSEIF ( simulated_time < 10000.0_wp )  THEN
286                WRITE( char_dum, '(I4)' )  INT( simulated_time )   
287             ELSEIF ( simulated_time < 100000.0_wp )  THEN
288                WRITE( char_dum, '(I5)' )  INT( simulated_time )   
289             ELSEIF ( simulated_time < 1000000.0_wp )  THEN
290                WRITE( char_dum, '(I6)' )  INT( simulated_time )
291             ELSEIF ( simulated_time < 10000000.0_wp )  THEN
292                WRITE( char_dum, '(I7)' )  INT( simulated_time )   
293             ELSEIF ( simulated_time < 100000000.0_wp )  THEN
294                WRITE( char_dum, '(I8)' )  INT( simulated_time )   
295             ELSEIF ( simulated_time < 1000000000.0_wp )  THEN
296                WRITE( char_dum, '(I9)' )  INT( simulated_time )
297             ENDIF             
298!
299!--          Copy HEADER file
300             CALL system('cp HEADER ' // TRIM( path ) // TRIM( char_dum ) //   & 
301                                      's_' // TRIM( variable_name ) // '.vtk')
302!--          Open VTK file.
303             OPEN ( file_id_out, FILE = TRIM( path ) // TRIM( char_dum ) //    & 
304                    's_' // TRIM( variable_name ) // '.vtk', FORM='FORMATTED', &
305                    POSITION = 'APPEND' )           
306          ENDIF
307!
308!--       Allocate and read array for variable data. 
309          ALLOCATE( var(1:ns(f)) )
310       
311          READ( file_id_in ) var
312!
313!--       Write variable data into output VTK file.
314          DO n = 1, ns(f)
315             WRITE( file_id_out, * ) var(n)
316          ENDDO
317!
318!--       Remember file position in binary file and close it.
319          filepos(f) = ftell( file_id_in )
320!           CALL FTELL( file_id_in, filepos(f) )
321       
322          CLOSE ( file_id_in )
323!
324!--       Deallocate temporary array for variable data.
325          DEALLOCATE( var )
326       
327       ENDDO
328!
329!--    After data for a variable for one time step is read, close the output
330!--    VTK file and go to next variable or timestep.
331       CLOSE ( file_id_out )
332!
333!--    If all files reached the end-of-file, exit the loop.
334       IF ( ALL( eof ) )  EXIT
335       
336    ENDDO
337!
338!-- Finally, remove HEADER file
339    CALL system( 'rm HEADER' )
340   
341 CONTAINS
342 
343!------------------------------------------------------------------------------!
344! Description:
345! ------------
346!> This subroutine read the namelist file.
347!------------------------------------------------------------------------------!
348    SUBROUTINE surface_output_parin
349       
350       IMPLICIT NONE
351       
352       INTEGER(iwp) ::  file_id_parin = 90
353       
354       NAMELIST /surface_output/  convert_average_data, cycle_number, num_pe,  &
355                                  path, run
356
357!
358!--    Open namelist file.
359       OPEN( file_id_parin, FILE='surface_output_parin',                       &
360             STATUS='OLD', FORM='FORMATTED')
361!
362!--    Read namelist.
363       READ ( file_id_parin, surface_output )
364!
365!--    Close namelist file.
366       CLOSE( file_id_parin )
367       
368    END SUBROUTINE surface_output_parin
369     
370!------------------------------------------------------------------------------!
371! Description:
372! ------------
373!> This subroutine creates the filename string of the treated binary file.
374!------------------------------------------------------------------------------!
375    SUBROUTINE surface_output_create_file_string
376       
377       IMPLICIT NONE
378       
379       CHARACTER(LEN=3) ::  char_av = ''
380       CHARACTER(LEN=4) ::  char_cycle = ''
381       
382!
383!--    Create substring for the cycle number.
384       IF ( cycle_number /= 0 )  THEN
385          IF ( cycle_number < 10 )  THEN
386             WRITE( char_cycle, '(I1)')  cycle_number
387             char_cycle = '.00' // TRIM( char_cycle )
388          ELSEIF ( cycle_number < 100 )  THEN
389             WRITE( char_cycle, '(I2)')  cycle_number
390             char_cycle = '.0' // TRIM( char_cycle )
391          ELSEIF ( cycle_number < 1000 )  THEN
392             WRITE( char_cycle, '(I3)')  cycle_number
393             char_cycle = '.' // TRIM( char_cycle )
394          ENDIF
395       ENDIF
396!
397!--    Create substring for averaged data.
398       IF ( convert_average_data )  char_av = '_av'
399!
400!--    Create substring for the processor id and combine all substrings.
401       IF ( f < 10 )  THEN
402          WRITE( char_dum, '(I1)')  f
403          myid_char = TRIM( char_av ) // '_surf_00000' // TRIM( char_dum ) //  &
404                      TRIM( char_cycle ) // file_suffix
405       ELSEIF ( f < 100     )  THEN
406          WRITE( char_dum, '(I2)')  f
407          myid_char = TRIM( char_av ) // '_surf_0000'  // TRIM( char_dum ) //  &
408                      TRIM( char_cycle ) // file_suffix
409       ELSEIF ( f < 1000    )  THEN
410          WRITE( char_dum, '(I3)')  f
411          myid_char = TRIM( char_av ) // '_surf_000'   // TRIM( char_dum ) //  &
412                      TRIM( char_cycle ) // file_suffix
413       ELSEIF ( f < 10000   )  THEN
414          WRITE( char_dum, '(I4)')  f
415          myid_char = TRIM( char_av ) // '_surf_00'    // TRIM( char_dum ) //  &
416                      TRIM( char_cycle ) // file_suffix
417       ELSEIF ( f < 100000  )  THEN
418          WRITE( char_dum, '(I5)')  f
419          myid_char = TRIM( char_av ) // '_surf_0'     // TRIM( char_dum ) //  &
420                      TRIM( char_cycle ) // file_suffix
421       ELSEIF ( f < 1000000 )  THEN
422          WRITE( char_dum, '(I6)')  f
423          myid_char = TRIM( char_av ) // '_surf_'      // TRIM( char_dum ) //  &
424                      TRIM( char_cycle ) // file_suffix
425       ENDIF
426       
427    END SUBROUTINE surface_output_create_file_string
428   
429 END PROGRAM surface_output_to_vtk
430
431
432
Note: See TracBrowser for help on using the repository browser.