Index: palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90
===================================================================
--- palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90	(revision 3805)
+++ palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90	(revision 3841)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add suffix to VTK output to indicate average data
+! 
+! 3805 2019-03-20 15:26:35Z raasch
 ! output format adjusted
 ! 
@@ -143,7 +146,11 @@
 !> Output is distinguished between instantaneous and time-averaged data. 
 !------------------------------------------------------------------------------!
- PROGRAM surface_output_to_vtk
+ PROGRAM surface_output_to_vtk_netcdf
  
     USE, INTRINSIC ::  ISO_C_BINDING
+    
+#if defined( __netcdf )  
+    USE NETCDF
+#endif
  
     USE posix_interface,                                                       &
@@ -171,5 +178,5 @@
     INTEGER, PARAMETER ::  OFFSET_KIND = C_SIZE_T !< unsigned integer for the C-interface
 
-    INTEGER(iwp) ::  cycle_number                 !< cycle number
+    INTEGER(iwp) ::  cycle_number = 0             !< cycle number
     INTEGER(iwp) ::  f                            !< running index over all binary files
     INTEGER(iwp) ::  file_id_in = 18              !< file unit for input binaray file    
@@ -189,4 +196,5 @@
     
     LOGICAL ::  convert_average_data = .FALSE.          !< namelist parameter to decide whether average or instantaneous data should be converted
+    LOGICAL ::  convert_to_netcdf    = .FALSE.          !< namelist parameter indicating that binary data is also converted to NetCDF files
     LOGICAL, DIMENSION(:), ALLOCATABLE      ::  eof     !< flag indicating that end of binary file is reached
     
@@ -217,5 +225,5 @@
 !-- VTK format.
     OPEN ( file_id_out_header, FILE = 'HEADER', STATUS = 'REPLACE',            &
-           FORM = 'FORMATTED' ) 
+           FORM = 'FORMATTED' )
 !
 !-- READ grid setup, i.e. the number and position of vertices and surface elements
@@ -223,6 +231,5 @@
 !-- 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. 
+!-- Hence, header information needs to be duplicated multiple times. 
 !-- Moreover, Paraview expects consecutive vertex and polygon data, which are 
 !-- all distributed over the binaray files. Hence, first read vertex data from 
@@ -245,9 +252,13 @@
        READ ( file_id_in )  ns(f)
        READ ( file_id_in )  ns_total
-
+!
+!--    If binary surface dater shall be converted into NetCDF, create the file.
+!--    Only for first call.
+       IF ( convert_to_netcdf  .AND.  f == 0 )  THEN
+          CALL surface_output_create_netcdf_file
+       ENDIF
 !
 !--    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.
@@ -404,11 +415,23 @@
              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' )           
+!--          Copy HEADER file and open VTK file
+             IF ( convert_average_data )  THEN
+                CALL system('cp HEADER ' //                                    &
+                            TRIM( path ) // TRIM( char_dum ) //                & 
+                            '_AV_' // 's_' // TRIM( variable_name ) // '.vtk' )
+                            
+                OPEN ( file_id_out, FILE = TRIM( path ) // TRIM( char_dum ) // & 
+                       '_AV_' // 's_' // TRIM( variable_name ) // '.vtk',      &
+                       FORM='FORMATTED', POSITION = 'APPEND' )  
+             ELSE
+                CALL system('cp HEADER ' //                                    &
+                            TRIM( path ) // TRIM( char_dum ) //                & 
+                            's_' // TRIM( variable_name ) // '.vtk' )
+                            
+                OPEN ( file_id_out, FILE = TRIM( path ) // TRIM( char_dum ) // & 
+                       's_' // TRIM( variable_name ) // '.vtk',                &
+                       FORM='FORMATTED', POSITION = 'APPEND' )  
+             ENDIF
+
           ENDIF
 !
@@ -459,6 +482,6 @@
        INTEGER(iwp) ::  file_id_parin = 90
        
-       NAMELIST /surface_output/  convert_average_data, cycle_number, num_pe,  &
-                                  path, run
+       NAMELIST /surface_output/  convert_average_data, convert_to_netcdf,     &
+                                  cycle_number, num_pe, path, run
 
 !
@@ -534,6 +557,103 @@
     END SUBROUTINE surface_output_create_file_string
     
- END PROGRAM surface_output_to_vtk
-
-
-
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine creates the NetCDF file and defines the dimensions
+!------------------------------------------------------------------------------!
+    SUBROUTINE surface_output_create_netcdf_file
+       
+       IMPLICIT NONE
+       
+       
+       CHARACTER(LEN=5)   ::  char_cycle = ''
+       CHARACTER(LEN=200) ::  nc_filename
+       
+!
+!--    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
+       ELSE
+          char_cycle = '.'
+       ENDIF
+#if defined( __netcdf )  
+!
+!-     Create file
+       nc_filename = TRIM( path ) // TRIM(run) // '_surf' //                   &
+                     TRIM( char_cycle ) // 'nc'
+                     
+       status_nc = NF90_CREATE( nc_filename, IOR( NF90_CLOBBER, NF90_NETCDF4 ),&
+                                nc_id )
+       CALL handle_error( 'create netcdf file' )
+!
+!--    Define dimensions
+       status_nc = NF90_DEF_DIM( nc_id, 'xs', ns_total, id_xs )
+       CALL handle_error( 'define dimension xs' )
+       status_nc = NF90_DEF_DIM( nc_id, 'ys', ns_total, id_ys )
+       CALL handle_error( 'define dimension ys' )
+       status_nc = NF90_DEF_DIM( nc_id, 'zs', ns_total, id_zs )
+       CALL handle_error( 'define dimension zs' )
+       status_nc = NF90_DEF_DIM( nc_id, 'azimuth', ns_total, id_az )
+       CALL handle_error( 'define dimension azimuth' )
+       status_nc = NF90_DEF_DIM( nc_id, 'zenith',  ns_total, id_ze )
+       CALL handle_error( 'define dimension zenith' )
+       status_nc = NF90_DEF_DIM( nc_id, 'time',  NF90_UNLIMITED, id_time )
+       CALL handle_error( 'define dimension time' )
+!
+!--    Define dimension variables
+       status_nc = NF90_DEF_VAR( nc_id, 'xs', NF90_FLOAT,                      &
+                                 (/ id_xs /), id_var_xs )
+       CALL handle_error( 'define variable xs' )
+       status_nc = NF90_DEF_VAR( nc_id, 'ys', NF90_FLOAT,                      &
+                                 (/ id_ys /), id_var_ys )
+       CALL handle_error( 'define variable ys' )
+       status_nc = NF90_DEF_VAR( nc_id, 'zs', NF90_FLOAT,                      &
+                                 (/ id_zs /), id_var_zs )
+       CALL handle_error( 'define variable zs' ) 
+       status_nc = NF90_DEF_VAR( nc_id, 'azimuth', NF90_FLOAT,                 &
+                                 (/ id_az /), id_var_az ) 
+       CALL handle_error( 'define variable azimuth' )         
+       status_nc = NF90_DEF_VAR( nc_id, 'zenith', NF90_FLOAT,                  &
+                                 (/ id_ze /), id_var_ze ) 
+       CALL handle_error( 'define variable zenith' )                          
+       status_nc = NF90_DEF_VAR( nc_id, 'time', NF90_FLOAT,                    &
+                                 (/ id_time /), id_var_time )
+       CALL handle_error( 'define variable time' )
+#endif        
+    END SUBROUTINE surface_output_create_netcdf_file
+    
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!>  Output of NetCDF error code
+!------------------------------------------------------------------------------!    
+    SUBROUTINE handle_error( action )
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*) ::  action
+       
+#if defined( __netcdf )       
+       IF ( status_nc /= NF90_NOERR )  THEN
+          PRINT*, TRIM( NF90_STRERROR( status_nc ) ) // ' -- ' // action 
+          STOP
+       ENDIF
+#endif
+
+    END SUBROUTINE handle_error
+    
+ END PROGRAM surface_output_to_vtk_netcdf
+
+
+
