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

Last change on this file since 3828 was 3805, checked in by raasch, 5 years ago

output format adjusted, unused variable removed

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