Index: palm/trunk/UTIL/Makefile_utilities
===================================================================
--- palm/trunk/UTIL/Makefile_utilities	(revision 3487)
+++ palm/trunk/UTIL/Makefile_utilities	(revision 3494)
@@ -23,4 +23,7 @@
 # -----------------
 # $Id$
+# Added make directive for surface-output post-processing
+# 
+# 3208 2018-08-27 13:10:50Z sward
 # Added make directive for agent_preprocessing tool
 # 
@@ -54,8 +57,10 @@
 PROG2 =  compare_palm_logs.x
 PROG3 =  agent_preprocessing
+PROG4 =  surface_output_to_vtk
 
 OBJS1 =  combine_plot_fields.o
 OBJS2 =  compare_palm_logs.o
 OBJS3 =  agent_preprocessing.o
+OBJS4 =  surface_output_to_vtk.o
 
 CC = cc
@@ -71,5 +76,5 @@
 
 
-all: $(PROG1) $(PROG2) $(PROG3)
+all: $(PROG1) $(PROG2) $(PROG3) $(PROG4)
 
 $(PROG1): $(OBJS1)
@@ -82,4 +87,7 @@
 	$(F90_SER) -o $(PROG3) $(OBJS3) $(LDFLAGS)
 
+$(PROG4): $(OBJS4)
+	$(F90_SER) -o $(PROG4) $(OBJS4) $(LDFLAGS)
+
 .f90.o:
 	$(F90_SER) $(F90FLAGS) $(COPT) -c $<
Index: palm/trunk/UTIL/surface_output_processing/surface_output_parin
===================================================================
--- palm/trunk/UTIL/surface_output_processing/surface_output_parin	(revision 3494)
+++ palm/trunk/UTIL/surface_output_processing/surface_output_parin	(revision 3494)
@@ -0,0 +1,18 @@
+&surface_output
+
+!
+!-- Number of processors of the run - equal the number of binary files with surface data.
+    num_pe = 4, 
+!
+!-- Path to the output binary files.
+    path = '~/palm/current_version/JOBS/test_urban/OUTPUT/', 
+!
+!-- Name of the run.
+    run = 'test_urban',
+!
+!-- Cycle number of the output data which should be processed.
+    cycle_number = 0, 
+!
+!-- Convert instantaneous or average data.
+    convert_average_data = .F., 
+/  
Index: palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90
===================================================================
--- palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90	(revision 3494)
+++ palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90	(revision 3494)
@@ -0,0 +1,424 @@
+!> @file surface_output_merge.f90
+!------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2018  Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Initial version
+! 
+! 3241 2018-09-12 15:02:00Z raasch
+!
+!
+! Authors:
+! --------
+! @author Matthias Suehring and Klaus Ketelsen
+!
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine combines surface output from PALM-subdomains into one file. 
+!> Output from every processor element is opened and read and output 
+!> from all processor elements are written into one file for each timestep 
+!> according to Paraview's VTK standard. 
+!> Output is distinguished between instantaneous and time-averaged data. 
+!------------------------------------------------------------------------------!
+ PROGRAM surface_output_merge
+
+    IMPLICIT NONE
+    
+    CHARACTER(LEN=4)   ::  char_time            !< string indicating simulated time
+    CHARACTER(LEN=4)   ::  file_suffix = '.bin' !< string which contain the suffix indicating surface data
+    
+    CHARACTER(LEN=10)  ::  char_dum             !< dummy string
+    
+    CHARACTER(LEN=30)  ::  myid_char            !< combined string indicating binary file
+    
+    CHARACTER(LEN=100) ::  path                 !< path to the binary data 
+    CHARACTER(LEN=100) ::  run                  !< name of the run
+    CHARACTER(LEN=100) ::  variable_name        !< name of the processed output variable    
+    
+    INTEGER(4)   ::  ftell                      !< intrinsic function, get current position in file
+    INTEGER(4)   ::  fseek                      !< intrinsic function, go to given position in file
+    INTEGER(4)   ::  ndum                       !< return parameter of intrinsic function fseek
+    
+    INTEGER, PARAMETER ::  iwp = 4              !< integer precision
+    INTEGER, PARAMETER ::  wp  = 8              !< float precision
+    
+    INTEGER(iwp) ::  cycle_number               !< cycle number
+    INTEGER(iwp) ::  f                          !< running index over all binary files
+    INTEGER(iwp) ::  file_id_in = 110           !< file unit for input binaray file    
+    INTEGER(iwp) ::  file_id_out = 20           !< file unit for output VTK file        
+    INTEGER(iwp) ::  file_id_out_header = 19    !< file unit for temporary header file
+    INTEGER(iwp) ::  length                     !< length of string on file
+    INTEGER(iwp) ::  n                          !< running index over all points and polygons
+    INTEGER(iwp) ::  npoints_total              !< total number of points
+    INTEGER(iwp) ::  ns_total                   !< total number of polygons
+    INTEGER(iwp) ::  num_pe                     !< number of processors of the run
+    
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  filepos !< current fileposition in binary file
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  npoints !< number of points/vertices in a binaray file
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns      !< number of surface elements in a binaray file 
+    
+    LOGICAL ::  convert_average_data = .FALSE.          !< namelist parameter to decide whether average or instantaneous data should be converted
+    LOGICAL, DIMENSION(:), ALLOCATABLE      ::  eof     !< flag indicating that end of binary file is reached
+    
+    REAL(wp)                              ::  simulated_time !< output time   
+    
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var            !< actual surface data
+    
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  points         !< point / vertex data
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  polygons       !< polygon data
+
+!
+!-- Read namelist.
+    CALL surface_output_merge_parin
+!
+!-- Allocate array which contains the file position in each output file,
+!-- in order to skip some reading.
+    ALLOCATE( eof( 0:num_pe-1 )     )
+    ALLOCATE( filepos( 0:num_pe-1 ) )
+    ALLOCATE( npoints( 0:num_pe-1 ) )
+    ALLOCATE( ns( 0:num_pe-1 )      )
+!
+!-- Initialize file position.
+    filepos = 0
+!
+!-- Open a temporary file which contains all necessary header information for the
+!-- VTK format.
+    OPEN ( file_id_out_header, FILE = 'HEADER', STATUS = 'REPLACE',            &
+           FORM = 'FORMATTED' ) 
+!
+!-- READ grid setup, i.e. the number and position of vertices and surface elements
+!-- and merge all this information into one file. Further, create all the required
+!-- header information of the VTK file.
+!-- Note, PARAVIEW expects one VTK file for each variable and each timestep. 
+!-- Hence, header information needs to be duplicated multiple times, which will be
+!-- be done later in a bash script. 
+!-- Moreover, Paraview expects consecutive vertex and polygon data, which are 
+!-- all distributed over the binaray files. Hence, first read vertex data from 
+!-- binary file, write this to the HEADER file, close the binary file and read
+!-- data from the next binary file and so on. This requires several openenings 
+!-- and closings of the binary files and temporarily storage of the 
+!-- file-positions.
+    DO  f = 0, num_pe - 1
+!
+!--    Create filename of the treated binary file.
+       CALL surface_output_merge_create_file_string
+!
+!--    Open file with surface output for processor f.
+       OPEN ( file_id_in + f, FILE = TRIM( path ) // TRIM( run ) //            &
+              TRIM( myid_char ), FORM='UNFORMATTED' )
+!
+!--    Read number of vertices / points and surface elements
+       READ ( file_id_in + f )  npoints(f)
+       READ ( file_id_in + f )  npoints_total
+       READ ( file_id_in + f )  ns(f)
+       READ ( file_id_in + f )  ns_total
+
+!
+!--    Allocate arrays where all the surfaces and vertices will be stored.
+       ALLOCATE( points(3,1:npoints(f))   )
+       
+!
+!--    Read polygon data and store them in a temporary file.
+       READ ( file_id_in + f )  points
+!
+!--    Obtain current file position. Will be stored for next file opening.
+       filepos(f) = ftell( file_id_in + f )
+!
+!--    Write header information. Only one time required.
+       IF ( f == 0 )  THEN
+          WRITE( file_id_out_header,'(A)' ) "# vtk DataFile Version 3.0"
+          WRITE( file_id_out_header,'(A,F8.2,A)' ) "legacy vtk File generated by PALM,  simulation time = xxx sec"
+          WRITE( file_id_out_header,'(A)') "ASCII"
+
+          WRITE( file_id_out_header,'(A)') "DATASET POLYDATA"
+          WRITE( file_id_out_header,'(A,I5,A)') "POINTS ", npoints_total, " float"
+       ENDIF      
+!
+!--    Write the vertex data into header file.
+       DO  n = 1, npoints(f)
+          WRITE( file_id_out_header, '(8F10.1)' )  points(1:3,n)
+       ENDDO
+!
+!--    Deallocate vertex data and close binary file.
+       DEALLOCATE( points )
+       
+       CLOSE ( file_id_in + f )
+    ENDDO
+!
+!-- Now, treat polygon data.
+    DO  f = 0, num_pe - 1
+!
+!--    Create filename of the treated binary file .   
+       CALL surface_output_merge_create_file_string
+!
+!--    Open file with surface output for processor f.
+       OPEN ( file_id_in + f, FILE = TRIM( path ) // TRIM( run ) //            &
+              TRIM( myid_char ), FORM='UNFORMATTED' )
+!
+!--    Move to last postion. 
+       ndum = fseek( file_id_in + f, filepos(f), 0 )
+!
+!--    Allocate array for polygon data
+       ALLOCATE( polygons(5,1:ns(f)) )
+!
+!--    Read polygon data and store them in a temporary file.
+       READ ( file_id_in + f )  polygons
+!
+!--    Obtain current file position after reading the local polygon data. 
+!--    Will be used for next file opening.
+       filepos(f) = ftell( file_id_in + f )
+!
+!--    Write further header information. Only one time required.
+       IF ( f == 0 )                                                           &
+          WRITE ( file_id_out_header, '(A,8I10)') "POLYGONS ",                 &
+                                                  ns_total, 5 * ns_total
+!
+!--    Write the polygons into the header file. Note, polygon data is of type
+!--    integer, as it just connects the point-indices which describe the given
+!--    surface element.
+       DO n = 1, ns(f)
+          WRITE ( file_id_out_header, '(8I10)' )  INT (polygons(1:5,n) )
+       ENDDO
+!
+!--    Deallocate array for polygon data and close the file.
+       DEALLOCATE( polygons )
+       
+       CLOSE ( file_id_in + f )
+       
+    ENDDO
+    
+    f = 0
+    CALL surface_output_merge_create_file_string
+!
+!-- Write further header information. Only once required.
+    WRITE ( file_id_out_header, '(A,I10)') "CELL_DATA ", ns_total
+    WRITE ( file_id_out_header, '(A,I10)') "SCALARS cell_scalars float 1 "
+    WRITE ( file_id_out_header, '(A,I10)') "LOOKUP_TABLE default "
+!
+!-- Header creation has now been completed. File can be closed.    
+    CLOSE( file_id_out_header )
+!
+!-- Now, read all the actual data from the surface output. Please note, Paraview
+!-- VTK format expects 1 file per variable and time step. 
+!-- The output file is only created once and includes the variable and the 
+!-- simulated time. 
+!-- In the binaray files, several variables and timesteps are stored, data for a 
+!-- given variable, however, is distributed over all binary files. Hence, read 
+!-- variable data for a given timestep from the binary file, write this data into
+!-- the target output file, remember file position in binary file and close it, 
+!-- open nex binary file and read the variable, and so on, until all variables
+!-- for all timesteps are processed. 
+    eof = .FALSE.
+    DO
+       DO  f = 0, num_pe - 1
+!
+!--       Clean up strings-
+          char_time = ''
+          variable_name = ''
+!
+!--       Create filename of the treated binary file.            
+          CALL surface_output_merge_create_file_string
+!
+!--       Open binary file with surface output for processor f.
+          OPEN ( file_id_in + f, FILE = TRIM( path ) // TRIM( run ) //         &
+                 TRIM( myid_char ), FORM='UNFORMATTED' )
+!
+!--       Move to last postion. 
+          ndum = fseek( file_id_in + f, filepos(f), 0 )
+!
+!--       Read string length and string indicating the output time step. 
+          READ ( file_id_in + f ) length
+          READ ( file_id_in + f ) char_time(1:length)
+!
+!--       If string for the output time indicates the end-of-file, set the eof 
+!--       flag and skip the read of the loop. 
+          IF ( char_time(1:length) == 'END' )  THEN
+             eof(f) = .TRUE. 
+             CLOSE ( file_id_in + f )
+             CYCLE
+          ENDIF
+!
+!--       Read output time, and variable name.
+          READ ( file_id_in + f ) simulated_time
+          READ ( file_id_in + f ) length
+          READ ( file_id_in + f ) variable_name(1:length)
+!
+!--       For first loop index, open the target output file. First create the
+!--       filename string. Further, copy HEADER file with the given filename 
+!--       string. The header information must be given in each VTK file!
+          IF ( f == 0 )  THEN
+             IF ( simulated_time < 10.0_wp )  THEN
+                WRITE( char_dum, '(I1)' )  INT( simulated_time ) 
+             ELSEIF ( simulated_time < 100.0_wp )  THEN
+                WRITE( char_dum, '(I2)' )  INT( simulated_time )
+             ELSEIF ( simulated_time < 1000.0_wp )  THEN
+                WRITE( char_dum, '(I3)' )  INT( simulated_time )
+             ELSEIF ( simulated_time < 10000.0_wp )  THEN
+                WRITE( char_dum, '(I4)' )  INT( simulated_time )   
+             ELSEIF ( simulated_time < 100000.0_wp )  THEN
+                WRITE( char_dum, '(I5)' )  INT( simulated_time )   
+             ELSEIF ( simulated_time < 1000000.0_wp )  THEN
+                WRITE( char_dum, '(I6)' )  INT( simulated_time )
+             ELSEIF ( simulated_time < 10000000.0_wp )  THEN
+                WRITE( char_dum, '(I7)' )  INT( simulated_time )   
+             ELSEIF ( simulated_time < 100000000.0_wp )  THEN
+                WRITE( char_dum, '(I8)' )  INT( simulated_time )   
+             ELSEIF ( simulated_time < 1000000000.0_wp )  THEN
+                WRITE( char_dum, '(I9)' )  INT( simulated_time )
+             ENDIF             
+!
+!--          Copy HEADER file
+             CALL system('cp HEADER ' // TRIM( path ) // TRIM( char_dum ) //   & 
+                                      's_' // TRIM( variable_name ) // '.vtk')
+!--          Open VTK file. 
+             OPEN ( file_id_out, FILE = TRIM( path ) // TRIM( char_dum ) //    & 
+                    's_' // TRIM( variable_name ) // '.vtk', FORM='FORMATTED', &
+                    POSITION='APPEND' )           
+          ENDIF
+!
+!--       Allocate and read array for variable data.  
+          ALLOCATE( var(1:ns(f)) )
+       
+          READ( file_id_in + f ) var
+!
+!--       Write variable data into output VTK file.
+          DO n = 1, ns(f)
+             WRITE( file_id_out, * ) var(n)
+          ENDDO
+!
+!--       Remember file position in binary file and close it.
+          filepos(f) = ftell( file_id_in + f )
+       
+          CLOSE ( file_id_in + f )
+!
+!--       Deallocate temporary array for variable data.
+          DEALLOCATE( var )
+       
+       ENDDO
+!
+!--    After data for a variable for one time step is read, close the output
+!--    VTK file and go to next variable or timestep.
+       CLOSE ( file_id_out )
+!
+!--    If all files reached the end-of-file, exit the loop. 
+       IF ( ALL( eof ) )  EXIT
+       
+    ENDDO
+!
+!-- Finally, remove HEADER file
+    CALL system( 'rm HEADER' )
+    
+ CONTAINS
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine read the namelist file.
+!------------------------------------------------------------------------------!
+    SUBROUTINE surface_output_merge_parin
+       
+       IMPLICIT NONE
+       
+       INTEGER(iwp) ::  file_id_parin = 90
+       
+       NAMELIST /surface_output/  convert_average_data, cycle_number, num_pe,  &
+                                  path, run
+
+!
+!--    Open namelist file.
+       OPEN( file_id_parin, FILE='surface_output_parin',                       &
+             STATUS='OLD', FORM='FORMATTED')
+!
+!--    Read namelist.
+       READ ( file_id_parin, surface_output )
+!
+!--    Close namelist file.
+       CLOSE( file_id_parin )
+       
+    END SUBROUTINE surface_output_merge_parin
+      
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine creates the filename string of the treated binary file.
+!------------------------------------------------------------------------------!
+    SUBROUTINE surface_output_merge_create_file_string
+       
+       IMPLICIT NONE
+       
+       CHARACTER(LEN=3) ::  char_av = ''
+       CHARACTER(LEN=4) ::  char_cycle = ''
+       
+!
+!--    Create substring for the cycle number.
+       IF ( cycle_number /= 0 )  THEN
+          IF ( cycle_number < 10 )  THEN
+             WRITE( char_cycle, '(I1)')  cycle_number
+             char_cycle = '.00' // TRIM( char_cycle )
+          ELSEIF ( cycle_number < 100 )  THEN
+             WRITE( char_cycle, '(I2)')  cycle_number
+             char_cycle = '.0' // TRIM( char_cycle )
+          ELSEIF ( cycle_number < 1000 )  THEN
+             WRITE( char_cycle, '(I3)')  cycle_number
+             char_cycle = '.' // TRIM( char_cycle )
+          ENDIF
+       ENDIF
+!
+!--    Create substring for averaged data. 
+       IF ( convert_average_data )  char_av = '_av'
+!
+!--    Create substring for the processor id and combine all substrings.
+       IF ( f < 10 )  THEN
+          WRITE( char_dum, '(I1)')  f
+          myid_char = TRIM( char_av ) // '_surf_00000' // TRIM( char_dum ) //  &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 100     )  THEN
+          WRITE( char_dum, '(I2)')  f
+          myid_char = TRIM( char_av ) // '_surf_0000'  // TRIM( char_dum ) //  &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 1000    )  THEN
+          WRITE( char_dum, '(I3)')  f
+          myid_char = TRIM( char_av ) // '_surf_000'   // TRIM( char_dum ) //  &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 10000   )  THEN
+          WRITE( char_dum, '(I4)')  f
+          myid_char = TRIM( char_av ) // '_surf_00'    // TRIM( char_dum ) //  &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 100000  )  THEN
+          WRITE( char_dum, '(I5)')  f
+          myid_char = TRIM( char_av ) // '_surf_0'     // TRIM( char_dum ) //  &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 1000000 )  THEN
+          WRITE( char_dum, '(I6)')  f
+          myid_char = TRIM( char_av ) // '_surf_'      // TRIM( char_dum ) //  &
+                      TRIM( char_cycle ) // file_suffix
+       ENDIF
+       
+    END SUBROUTINE surface_output_merge_create_file_string
+    
+ END PROGRAM surface_output_merge
+
+
+
