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

Last change on this file since 4338 was 3842, checked in by suehring, 6 years ago

remove redundant VTK netcdf option in postprocessor

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