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

!> @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: surface_output_to_vtk.f90 3841 2019-03-29 11:38:44Z suehring $
! Add suffix to VTK output to indicate average data
! 
! 3805 2019-03-20 15:26:35Z raasch
! output format adjusted
! 
! 3755 2019-02-19 17:03:19Z suehring
! - Change format description to avoid *** for larger domains
! - switch back from to non-standard Fortran ftell and fseek since posix
!   interface do not work any more for some unknown reason
! 
! 3523 2018-11-13 16:09:31Z suehring
! Implement interface for posix conform C-systemcalls in order to replace
! non-standard FORTRAN functions ftell and fseek. 
! 
! 3496 2018-11-06 15:59:50Z suehring
! Use subroutine call for fseek instead of function call. gfortran has some
! problems with this.
! 
! 3494 2018-11-06 14:51:27Z suehring
! Initial version
!
! Authors:
! --------
! @author Matthias Suehring and Klaus Ketelsen
!
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Interface to posix systemcalls. The module allows to call the following
!> C system calls from FORTRAN: ftell, lseek.
!------------------------------------------------------------------------------!
 MODULE posix_interface
 
    INTERFACE
       INTEGER (C_SIZE_T) FUNCTION C_LSEEK (fd, offset, whence)                &
                                   BIND(C, NAME='lseek')
          
          USE ISO_C_BINDING
          
          IMPLICIT NONE
          
          INTEGER(KIND=C_INT),    VALUE ::  fd
          INTEGER(KIND=C_SIZE_T), VALUE ::  offset
          INTEGER(KIND=C_INT),    VALUE ::  whence
          
       END FUNCTION C_LSEEK
       
    END INTERFACE
    
    INTERFACE
       INTEGER (C_SIZE_T) FUNCTION C_FTELL ( fd )                              &
                                   BIND(C, NAME='ftell')
          
          USE ISO_C_BINDING
          
          IMPLICIT NONE
          
          INTEGER(KIND=C_INT),    VALUE ::  fd

       END FUNCTION C_FTELL
       
    END INTERFACE
    
    PUBLIC posix_lseek, posix_ftell
    
    CONTAINS
!
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Interface for the C-routine lseek. 
!------------------------------------------------------------------------------! 
    SUBROUTINE posix_lseek( fid, offset )
    
       USE ISO_C_BINDING
    
       IMPLICIT NONE

       INTEGER,INTENT(IN)                     ::  fid    !< file unit
       INTEGER(KIND=C_SIZE_T),INTENT(IN)      ::  offset !< file offset from the beginning

       INTEGER(KIND=C_INT)                    ::  my_fid !< file unit
       INTEGER(KIND=C_SIZE_T)                 ::  retval !< return value
       INTEGER(KIND=C_INT)                    ::  whence !< mode, here start from the beginning

       my_fid = fid
       whence = 0

       retval = C_LSEEK( fid, offset, whence )
       
    END SUBROUTINE posix_lseek
!
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Interface for the C-routine ftell. 
!------------------------------------------------------------------------------!
    SUBROUTINE posix_ftell( fid, filepos )
    
       USE ISO_C_BINDING
    
       IMPLICIT NONE

       INTEGER,INTENT(IN)      ::  fid                     !< file unit
       INTEGER(KIND=C_SIZE_T), INTENT(INOUT)   ::  filepos !< file position

       filepos = C_FTELL( fid )
       
    END SUBROUTINE posix_ftell

 END MODULE posix_interface
!
!------------------------------------------------------------------------------!
! 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_to_vtk_netcdf
 
    USE, INTRINSIC ::  ISO_C_BINDING
    
#if defined( __netcdf )  
    USE NETCDF
#endif
 
    USE posix_interface,                                                       &
        ONLY:  posix_ftell, posix_lseek

    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)   ::  ndum                       !< return parameter of intrinsic function fseek   
       
    INTEGER, PARAMETER ::  iwp = 4                !< integer precision
    INTEGER, PARAMETER ::  wp  = 8                !< float precision
    INTEGER, PARAMETER ::  OFFSET_KIND = C_SIZE_T !< unsigned integer for the C-interface

    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    
    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(OFFSET_KIND),DIMENSION(:), ALLOCATABLE ::  filepos !< current fileposition in binary file
    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 ::  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
    
    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
    
    logical :: flag

!
!-- Read namelist.
    CALL surface_output_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. 
!-- 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_create_file_string
!
!--    Open file with surface output for processor f.
       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
              TRIM( myid_char ), FORM = 'UNFORMATTED' )
!
!--    Read number of vertices / points and surface elements
       READ ( file_id_in )  npoints(f)
       READ ( file_id_in )  npoints_total
       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.
       READ ( file_id_in )  points
!
!--    Obtain current file position. Will be stored for next file opening.
!        CALL posix_ftell( file_id_in, filepos(f) )
       filepos(f) = ftell( file_id_in )
!
!--    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,I12,A)') "POINTS ", npoints_total, " float"
       ENDIF      
!
!--    Write the vertex data into header file.
       DO  n = 1, npoints(f)
          WRITE( file_id_out_header, '(8F15.1)' )  points(1:3,n)
       ENDDO
!
!--    Deallocate vertex data and close binary file.
       DEALLOCATE( points )
       
       CLOSE ( file_id_in )
    ENDDO
!
!-- Now, treat polygon data.
    DO  f = 0, num_pe - 1
!
!--    Create filename of the treated binary file .   
       CALL surface_output_create_file_string
!
!--    Open file with surface output for processor f.
       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
              TRIM( myid_char ), FORM = 'UNFORMATTED' )
!
!--    Move to last postion. 
!        CALL posix_lseek( file_id_in, filepos(f) )
       CALL FSEEK( file_id_in, filepos(f), 0, ndum )
!
!--    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 )  polygons
!
!--    Obtain current file position after reading the local polygon data. 
!--    Will be used for next file opening.
!        CALL posix_ftell( file_id_in, filepos(f) )
       filepos(f) = ftell( file_id_in )
!
!--    Write further header information. Only one time required.
       IF ( f == 0 )                                                           &
          WRITE ( file_id_out_header, '(A,2I10)') "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, '(5I18)' )  INT( polygons(1:5,n) )
       ENDDO
!
!--    Deallocate array for polygon data and close the file.
       DEALLOCATE( polygons )
       
       CLOSE ( file_id_in )
       
    ENDDO
    
    f = 0
    CALL surface_output_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_create_file_string
!
!--       Open binary file with surface output for processor f.
          OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //         &
                 TRIM( myid_char ), FORM = 'UNFORMATTED' )
!
!--       Move to last postion. 
!           CALL posix_lseek( file_id_in, filepos(f) )
          CALL FSEEK( file_id_in, filepos(f), 0, ndum )
!
!--       Read string length and string indicating the output time step. 
          READ ( file_id_in ) length
          READ ( file_id_in ) 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 )
             CYCLE
          ENDIF
!
!--       Read output time, and variable name.
          READ ( file_id_in ) simulated_time
          READ ( file_id_in ) length
          READ ( file_id_in ) 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 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
!
!--       Allocate and read array for variable data.  
          ALLOCATE( var(1:ns(f)) )
       
          READ( file_id_in ) 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.
!           CALL posix_ftell( file_id_in, filepos(f) )
          filepos(f) = ftell( file_id_in )
       
          CLOSE ( file_id_in )
!
!--       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_parin
       
       IMPLICIT NONE
       
       INTEGER(iwp) ::  file_id_parin = 90
       
       NAMELIST /surface_output/  convert_average_data, convert_to_netcdf,     &
                                  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_parin
      
!------------------------------------------------------------------------------!
! Description:
! ------------
!> This subroutine creates the filename string of the treated binary file.
!------------------------------------------------------------------------------!
    SUBROUTINE surface_output_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_create_file_string
    
    
!------------------------------------------------------------------------------!
! 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