Index: palm/trunk/SCRIPTS/.palm.iofiles
===================================================================
--- palm/trunk/SCRIPTS/.palm.iofiles	(revision 3726)
+++ palm/trunk/SCRIPTS/.palm.iofiles	(revision 3727)
@@ -59,4 +59,6 @@
 SURFACE_DATA_BIN*          out:tr   *         $base_data/$run_identifier/OUTPUT         _surf      bin
 SURFACE_DATA_AV_BIN*       out:tr   *         $base_data/$run_identifier/OUTPUT         _av_surf   bin
+SURFACE_DATA_NETCDF*       out:tr   *         $base_data/$run_identifier/OUTPUT         _surf      nc
+SURFACE_DATA_AV_NETCDF*    out:tr   *         $base_data/$run_identifier/OUTPUT         _av_surf   nc
 #
 VIRTUAL_MEAS_BIN*          out   *            $base_data/$run_identifier/OUTPUT         _vmeas     bin
Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3726)
+++ palm/trunk/SOURCE/Makefile	(revision 3727)
@@ -25,15 +25,18 @@
 # -----------------
 # $Id$
+# surface_data_output_mod depends on netcdf_interface_mod
+#
+# 3700 2019-01-26 17:03:42Z knoop
 # Rename surface_output_mod into surface_data_output_mod
-# 
+#
 # 3637 2018-12-20 01:51:36Z knoop
 # Implementation of the PALM module interface
-# 
+#
 # 3634 2018-12-18 12:31:28Z knoop
 # OpenACC port for SPEC
-# 
+#
 # 3579 2018-11-29 15:32:39Z suehring
 # Dependency for check_parameters on nesting_offl_mod added
-# 
+#
 # 3569 2018-11-27 17:03:40Z kanani
 # dom_dwd_user, Schrempf:
@@ -43,22 +46,22 @@
 # 3525 2018-11-14 16:06:14Z kanani
 # Changes related to clean-up of biometeorology (dom_dwd_user)
-# 
+#
 # 3522 2018-11-13 12:14:36Z suehring
 # Dependencies for virtual measurement module added
-# 
+#
 # 3494 2018-11-06 14:51:27Z suehring
 # Surface output revised
-# 
+#
 # 3474 2018-10-30 21:07:39Z kanani
 # Add virtual measurement module
-# 
+#
 # 3472 2018-10-30 20:43:50Z suehring
 # Add indoor model (kanani, srissman, tlang),
 # minor formatting
-# 
+#
 # 3467 2018-10-30 19:05:21Z suehring
 # Implementation of a new aerosol module salsa.
-# 
-# 
+#
+#
 # 3458 2018-10-30 14:51:23Z kanani
 # from chemistry branch r3443, banzhafs, Russo, forkel, basit:
@@ -67,23 +70,23 @@
 # Added chemistry emission module
 # chemistry_model_mod added to flow_statistics
-# 
+#
 # 3448 2018-10-29 18:14:31Z kanani
 # Adjustment of biometeorology dependencies
-# 
+#
 # 3436 2018-10-26 18:35:15Z gronemeier
 # Add surface_mod to user_data_output_mask
-# 
+#
 # 3435 2018-10-26 18:25:44Z gronemeier
 # - Add surface_mod to data_output_mask
 # - Add chemistry_model_mod and surface_mod to init_masks
-# 
+#
 # 3421 2018-10-24 18:39:32Z gronemeier
 # Add netcdf_data_input_mod to netcdf_interface_mod
 # bugfix: add dependencies to chemistry_model_mod
 # Add module for surface data output
-# 
+#
 # 3381 2018-10-19 13:09:06Z raasch
 # dependencies for ocean_mod fixed
-# 
+#
 # 3355 2018-10-16 14:03:34Z knoop
 # Add module for offline nesting;
@@ -91,13 +94,13 @@
 # Bugfix, missing dependency for turbulence generator in init_3d_model;
 # Some formatting ajdustments
-# 
+#
 # 3343 2018-10-15 10:38:52Z suehring
 # (from branch resler)
 # Add biometeorology
-# 
-# 
+#
+#
 # 3322 2018-10-09 10:02:39Z kanani
 # Formatting and cleanup
-# 
+#
 # 3298 2018-10-02 12:21:11Z kanani
 # Added missing dependencies and replaced blanks with tabs (forkel)
@@ -108,68 +111,68 @@
 # changes related to modularization of the ocean mode,
 # bugfix: dependency to advec_ws was missed in chemistry_model_mod
-# 
+#
 # 3274 2018-09-24 15:42:55Z knoop
 # Added palm dependency of multi_agent_system_mod, because of mas_last_actions
 # call at the end of palm run
-# 
+#
 # 3167 2018-07-24 18:17:30Z suehring
 # Bugfix, add missing dependencies for multi-agent system
-# 
+#
 # 3159 2018-07-20 11:20:01Z sward
 # Added multi agent system
-# 
+#
 # 3130 2018-07-16 11:08:55Z gronemeier
 # add surface_layer_fluxes_mod to turbulence_closure_mod
-# 
+#
 # 3129 2018-07-16 07:45:13Z gronemeier
 # add turbulence_closure_mod to parin
-# 
+#
 # 2963 2018-04-12 14:47:44Z suehring
-# Introduce index for vegetation/wall, pavement/green-wall and water/window 
+# Introduce index for vegetation/wall, pavement/green-wall and water/window
 # surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
-# 
+#
 # 2955 2018-04-09 15:14:01Z suehring
 # Add log-points to measure CPU time of NetCDF data input.
-# 
+#
 # 2938 2018-03-27 15:52:42Z suehring
 # No initialization of child domains via dynamic input file, except for soil
 # moisture and temperature
 # Apply turbulence generator at non-cyclic lateral boundary in nesting case
-# 
+#
 # 2936 2018-03-27 14:49:27Z suehring
 # Added dependencies for parent and child synchronization
-# 
+#
 # 2921 2018-03-22 15:05:23Z Giersch
 # date_and_time_mod dependency has been added to read/write_restart_data_mod
-# 
+#
 # 2918 2018-03-21 15:52:14Z gronemeier
-# read/write_3d_binary and read/write_var_list has been removed, 
+# read/write_3d_binary and read/write_var_list has been removed,
 # read/write_restart_data_mod, wrd_write_string and
-# user_read/write_restart_data_mod has been added, dependencies with respect to 
+# user_read/write_restart_data_mod has been added, dependencies with respect to
 # the aforementioned routines have been added/removed
 #
 # 2847 2018-03-02 21:45:58Z suehring
 # Changed format and enforced sorting
-# 
+#
 # 2817 2018-02-19 16:32:21Z knoop
 # Preliminary gust module interface implemented
-# 
+#
 # 2802 2018-02-14 16:21:39Z thiele
 # Changed lpm from subroutine to module.
 # Introduce particle transfer in nested models.
-# 
+#
 # 2773 2018-01-30 14:12:54Z suehring
 # Nesting of chemical species
-# 
+#
 # 2718 2018-01-02 08:49:38Z maronga
 # Corrected "Former revisions" section
-# 
+#
 # 2697 2017-12-14 17:57:20Z kanani
 # Bugfix, missing dependencies
-# 
+#
 # 2696 2017-12-14 17:12:51Z kanani
 # Change in file header (GPL part)
 # Implementation of uv exposure model (FK)
-# Bugfix, removed loop dependcy for vertical_nesting_mod and 
+# Bugfix, removed loop dependcy for vertical_nesting_mod and
 # turbulence_closure_mod, added depencies for vertical_nesting_mod (TG)
 # implemented turbulence_closure_mod (TG)
@@ -179,13 +182,13 @@
 # For LSM, add dependency on calc_mean_profile (??)
 # poismg_noopt modularized and renamed into poismg_noopt_mod
-# add dependencies for netcdf_data_input_mod, calc_mean_profile, 
+# add dependencies for netcdf_data_input_mod, calc_mean_profile,
 # radiation_module_mod, land_surface_model_mod (MS)
-# 
+#
 # 2608 2017-11-13 14:04:26Z schwenkel
-# Added diagnostic_quantities_mod 
-# 
+# Added diagnostic_quantities_mod
+#
 # 2600 2017-11-01 14:11:20Z raasch
 # comment line concerning bound checks removed
-# 
+#
 # 2599 2017-11-01 13:18:45Z hellstea
 # virtual_flight_mod, synthetic_turbulence_generator_mod and
@@ -194,40 +197,40 @@
 # 2563 2017-10-19 15:36:10Z Giersch
 # wind_turbine_model_mod and synthetic_turbulence_generator_mod were added to
-# write_var_list and virtual_flight_mod was deleted from read_var_list 
+# write_var_list and virtual_flight_mod was deleted from read_var_list
 #
 # 2544 2017-10-13 18:09:32Z maronga
-# Added date_and_time_mod 
-# 
+# Added date_and_time_mod
+#
 # 2371 2017-08-24 13:01:17Z kanani
 # Corrected dependencies for vertical_nesting_mod
-# 
+#
 # 2370 2017-08-23 06:11:43Z raasch
 # dependency bugfix for synthetic_turbulence_generator
-# 
+#
 # 2365 2017-08-21 14:59:59Z kanani
 # Added dependencies for vertical_nesting_mod
-# 
+#
 # 2339 2017-08-07 13:55:26Z gronemeier
 # corrected timestamp in header
-# 
+#
 # 2338 2017-08-07 12:15:38Z gronemeier
 # Modularize 1D model
-# 
+#
 # 2320 2017-07-21 12:47:43Z suehring
 # -ls_forcing nudging
 # +large_scale_forcing_nudging
-# 
+#
 # 2318 2017-07-20 17:27:44Z suehring
 # Add further dependencies on surface_mod
-# 
+#
 # 2317 2017-07-20 17:27:19Z suehring
 # Added time_integration_spinup
-# 
+#
 # 2269 2017-06-09 11:57:32Z suehring
 # Add dependency in read_3d_binary
-# 
+#
 # 2263 2017-06-08 14:59:01Z schwenkel
 # Implemented splitting and merging algorithm
-# 
+#
 # 2259 2017-06-08 09:09:11Z gronemeier
 # Implemented synthetic turbulence generator
@@ -235,12 +238,12 @@
 # 2256 2017-06-07 13:58:08Z suehring
 # Remove ring dependency in init_pegrid
-# 
+#
 # 2238 2017-05-31 16:49:16Z suehring
 # Bugfix, further missing dependency on surface_mod
-# 
+#
 # 2237 2017-05-31 10:34:53Z suehring
-# Bugfix, add dependencies on surface_mod for surface_coupler, 
-# plant_canopy_model_mod and ls_forcing_mod 
-# 
+# Bugfix, add dependencies on surface_mod for surface_coupler,
+# plant_canopy_model_mod and ls_forcing_mod
+#
 # 2233 2017-05-30 18:08:54Z suehring
 #
@@ -248,5 +251,5 @@
 # +dependencies for surface_mod
 # -wall_fluxes
-# 
+#
 # 2130 2017-01-24 16:25:39Z raasch
 # dependency for timestep updated
@@ -254,13 +257,13 @@
 # 2118 2017-01-17 16:38:49Z raasch
 # -cuda_fft_interfaces_mod
-# 
+#
 # 2050 2016-11-08 15:00:55Z gronemeier
 # Implement turbulent outflow method
-# 
+#
 # 2007 2016-08-24 15:47:17Z kanani
 # urban surface module added,
 # cleaned up some lines (compiler flags/options), which were accidentally
 # added in rev1938
-# 
+#
 # 1998 2016-08-20 18:45:34Z knoop
 # Bugfix: added netcdf_interface to dependency list for user_init_land_surface
@@ -268,8 +271,8 @@
 # 1986 2016-08-10 14:07:17Z gronemeier
 # POSIX-calls module added
-# 
+#
 # 1972 2016-07-26 07:52:02Z maronga
 # Removed some dependencies due to further modularization of land surface model
-# 
+#
 # 1957 2016-07-07 10:43:48Z suehring
 # flight module added
@@ -279,9 +282,9 @@
 #
 # 1938 2016-06-13 15:26:05Z hellstea
-# Some dependency errors corrected 
-# 
+# Some dependency errors corrected
+#
 # 1934 2016-06-13 09:46:57Z hellstea
 # poismg renamed poismg_noopt, poismg_fast_mod renamed poismg_mod
-# 
+#
 # 1914 2016-05-26 14:44:07Z witha
 # Added wind_turbine_model_mod
@@ -294,5 +297,5 @@
 #
 # 1850 2016-04-08 13:29:27Z maronga
-# Adapted for modularization of microphysics 
+# Adapted for modularization of microphysics
 # Several files renamed --> _mod
 # Bugfix for previous commit
@@ -309,16 +312,16 @@
 # 1833 2016-04-07 14:23:03Z raasch
 # spectrum renamed spectra_mod, depencies for spectra changed
-# 
+#
 # 1826 2016-04-07 12:01:39Z maronga
 # Renamed radiation_model to radiation_model_mod.
 # Renamed plant_canopy_model to plant_canopy_model_mod.
-# 
+#
 # 1822 2016-04-07 07:49:42Z hoffmann
 # Tails removed. lpm_release_set removed. calc_precipitation, impact_of_latent_heat
-# removed. 
+# removed.
 #
 # 1817 2016-04-06 15:44:20Z maronga
 # Renamed land_surface_model to land_surface_model_mod.
-# 
+#
 # 1808 2016-04-05 19:44:00Z raasch
 # -local_flush, -local_getenv
@@ -346,5 +349,5 @@
 #
 # 1762 2016-02-25 12:31:13Z hellstea
-# +pmc_interface, +pmc routines 
+# +pmc_interface, +pmc routines
 #
 # 1747 2016-02-08 12:25:53Z raasch
@@ -352,10 +355,10 @@
 #
 # 1691 2015-10-26 16:17:44Z maronga
-# Replaced prandtl_fluxes with surface_layer_fluxes. Added radiation_model to 
+# Replaced prandtl_fluxes with surface_layer_fluxes. Added radiation_model to
 # prognostic_equations
-# 
+#
 # 1585 2015-04-30 07:05:52Z maronga
 # Added user_init_radiation.f90
-# 
+#
 # 1575 2015-03-27 09:56:27Z raasch
 # +poismg_fast
@@ -363,5 +366,5 @@
 # 1551 2015-03-03 14:18:16Z maronga
 # Bugfix: further adjustments for the land surface model and radiation model
-# 
+#
 # 1517 2015-01-07 19:12:25Z hoffmann
 # advec_s_bc added to prognostic_equations
@@ -369,16 +372,16 @@
 # 1500 2014-12-03 17:42:41Z maronga
 # Bugfix: missing adjustments for land surface model and radiation model
-# 
+#
 # 1496 2014-12-02 17:25:50Z maronga
-# Added land surface model and radiation model files: land_surface_model, 
+# Added land surface model and radiation model files: land_surface_model,
 # radiation_model, user_init_land_surface
-# 
+#
 # 1484 2014-10-21 10:53:05Z kanani
 # plant_canopy_model-dependency added for check_parameters, header, init_3d_model,
 # package_parin, read_var_list, user_init_plant_canopy, write_var_list
-# 
+#
 # 1444 2014-08-02 20:10:32Z letzel
 # bugfix: cpulog added to lpm_advec
-# 
+#
 # 1404 2014-05-14 09:01:39Z keck
 # bugfix: dependencies added for progress_bar
@@ -386,11 +389,11 @@
 # 1402 2014-05-09 14:25:13Z raasch
 # progress_bar added
-# 
+#
 # 1400 2014-05-09 14:03:54Z knoop
 # Added new module random_generator_parallel
-# 
+#
 # 1380 2014-04-28 12:40:45Z heinze
 # bugfix: mod_particle_attributes  added to check_open
-# nudging added to time_integration 
+# nudging added to time_integration
 #
 # 1374 2014-04-25 12:55:07Z raasch
@@ -400,5 +403,5 @@
 # Added new module calc_mean_profile, previously in module buoyancy,
 # removed buoyancy dependency from nudging
-# 
+#
 # 1363 2014-04-17 12:28:49Z keck
 # bugfix: cpulog added to lpm_pack_arrays
@@ -406,7 +409,7 @@
 # 1361 2014-04-16 15:17:48Z hoffmann
 # cpulog added to microphysics
-# 
+#
 # 1359 2014-04-11 17:15:14Z hoffmann
-# mod_particle_attributes added, lpm_sort_arrays removed, 
+# mod_particle_attributes added, lpm_sort_arrays removed,
 # lpm_extend_particle_array removed
 #
@@ -715,5 +718,5 @@
 
 CC = cc
-CFLAGS = -O 
+CFLAGS = -O
 
 F90 =
@@ -1628,4 +1631,6 @@
 	mod_kinds.o \
 	modules.o \
+	netcdf_data_input_mod.o \
+	netcdf_interface_mod.o \
 	surface_mod.o
 swap_timelevel.o: \
Index: palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3726)
+++ palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3727)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! make several routines publicly available
+! 
+! 3701 2019-01-26 18:57:21Z knoop
 ! Statement added to prevent compiler warning about unused variable
 ! 
@@ -599,4 +602,8 @@
     END INTERFACE netcdf_create_file
 
+    INTERFACE netcdf_create_global_atts
+       MODULE PROCEDURE netcdf_create_global_atts
+    END INTERFACE netcdf_create_global_atts
+
     INTERFACE netcdf_create_var
        MODULE PROCEDURE netcdf_create_var
@@ -619,5 +626,6 @@
     END INTERFACE netcdf_open_write_file
 
-    PUBLIC netcdf_create_file, netcdf_define_header,                           &
+    PUBLIC netcdf_create_att, netcdf_create_dim, netcdf_create_file,           &
+           netcdf_create_global_atts, netcdf_create_var, netcdf_define_header, &
            netcdf_handle_error, netcdf_open_write_file
 
@@ -6756,5 +6764,5 @@
           CALL message( 'define_netcdf_header', 'PA0259', 0, 0, 0, 6, 0 )
 
-          
+
        CASE DEFAULT
 
Index: palm/trunk/SOURCE/surface_data_output_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_data_output_mod.f90	(revision 3726)
+++ palm/trunk/SOURCE/surface_data_output_mod.f90	(revision 3727)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enable NetCDF output for surface data (suehring, gronemeier)
+! 
+! 3691 2019-01-23 09:57:04Z suehring
 ! Add output of surface-parallel flow speed 
 ! 
@@ -77,5 +80,7 @@
        
    USE control_parameters,                                                     &
-       ONLY:  surface_output
+       ONLY:  coupling_char, data_output_during_spinup, end_time,              &
+              message_string, run_description_header, simulated_time_at_begin, &
+              spinup_time, surface_output
        
    USE grid_variables,                                                         &
@@ -85,4 +90,15 @@
        ONLY: nxl, nxr, nys, nyn, nzb, nzt
        
+#if defined( __netcdf )
+   USE NETCDF
+#endif
+
+   USE netcdf_data_input_mod,                                                  &
+       ONLY:  init_model
+
+   USE netcdf_interface,                                                       &
+       ONLY:  netcdf_create_att, netcdf_create_dim, netcdf_create_file,        &
+              netcdf_create_global_atts, netcdf_create_var, netcdf_handle_error
+       
    USE pegrid
        
@@ -100,6 +116,15 @@
       INTEGER(iwp) ::  npoints_total  !< total number of points / vertices which define a surface element
       
+      INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  s        !< coordinate for NetCDF output, number of the surface element
+      
       REAL(wp) ::  fillvalue = -9999.9_wp !< fillvalue for surface elements which are not defined
       
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  azimuth  !< azimuth orientation coordinate for NetCDF output 
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  es_utm   !< E-UTM coordinate for NetCDF output 
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  ns_utm   !< E-UTM coordinate for NetCDF output 
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  xs       !< x-coordinate for NetCDF output 
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  ys       !< y-coordinate for NetCDF output 
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zs       !< z-coordinate for NetCDF output 
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zenith   !< zenith orientation coordinate for NetCDF output 
       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_out  !< output variables
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_av   !< variables used for averaging
@@ -109,11 +134,33 @@
    
    CHARACTER(LEN=100), DIMENSION(300)       ::  data_output_surf = ' '  !< namelist variable which describes the output variables
-   CHARACTER(LEN=100), DIMENSION(0:1,300)   ::  dosurf = ' '            !< internal variable which describes the output variables and separates 
-                                                                        !< averaged from non-averaged output
+   CHARACTER(LEN=100), DIMENSION(0:1,300)   ::  dosurf = ' '            !< internal variable which describes the output variables 
+                                                                        !<  and separates averaged from non-averaged output
+   CHARACTER(LEN=100), DIMENSION(0:1,300)   ::  dosurf_unit = ' '       !< internal variable which holds the unit of the given output variable
    
-   INTEGER(iwp) ::  average_count_surf = 0                !< number of ensemble members used for averaging
-   INTEGER(iwp) ::  dosurf_no(0:1) = 0                    !< number of surface output quantities
+   INTEGER(iwp) ::  average_count_surf = 0   !< number of ensemble members used for averaging
+   INTEGER(iwp) ::  dosurf_no(0:1) = 0       !< number of surface output quantities
+
+   INTEGER(iwp) :: nc_stat                   !< error code for netcdf routines
+   INTEGER(iwp), SAVE ::  oldmode            !< save old set-fill-mode of netcdf file (not needed, but required for routine call)
+
+   INTEGER(iwp), DIMENSION(0:1) ::  id_dim_s_surf         !< netcdf ID for dimension s
+   INTEGER(iwp), DIMENSION(0:1) ::  id_dim_time_surf      !< netcdf ID for dimension time
+   INTEGER(iwp), DIMENSION(0:1) ::  id_set_surf           !< netcdf ID for file
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_etum_surf      !< netcdf ID for variable Es_UTM
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_nutm_surf      !< netcdf ID for variable Ns_UTM
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_time_surf      !< netcdf ID for variable time
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_s_surf         !< netcdf ID for variable s
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_xs_surf        !< netcdf ID for variable xs
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_ys_surf        !< netcdf ID for variable ys
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_zs_surf        !< netcdf ID for variable zs
+   INTEGER(iwp), DIMENSION(0:1) ::  dosurf_time_count = 0 !< count of output time steps
+   INTEGER(iwp), DIMENSION(0:1) ::  ntdim_surf            !< number of output time steps
+   
+   INTEGER(iwp), DIMENSION(0:1,300) ::  id_var_dosurf     !< netcdf ID for output variables
 
    LOGICAL :: first_output(0:1) = .FALSE.                 !< true if first output was already called
+   LOGICAL :: to_netcdf = .TRUE.                          !< flag indicating parallel NetCDF output 
+   LOGICAL :: to_vtk = .TRUE.                             !< flag indicating binary surface-data output that can be further 
+                                                          !< processed to VTK format
 
    REAL(wp) ::  averaging_interval_surf  = 9999999.9_wp   !< averaging interval
@@ -125,5 +172,4 @@
    REAL(wp) ::  time_dosurf_av = 0.0_wp                   !< internal counter variable to check for averaged data output
 
-   
    TYPE(surf_out) ::  surfaces      !< variable which contains all required output information
    
@@ -178,4 +224,9 @@
       IMPLICIT NONE
 
+      CHARACTER (LEN=100)  :: filename            !< name of output file
+      CHARACTER (LEN=80)   ::  time_average_text  !< string written to file attribute time_avg
+      CHARACTER (LEN=4000) ::  var_list           !< list of variables written to NetCDF file
+
+      INTEGER(iwp) ::  av                !< flag for averaged (=1) and non-averaged (=0) data
       INTEGER(iwp) ::  i                 !< grid index in x-direction, also running variable for counting non-average data output
       INTEGER(iwp) ::  j                 !< grid index in y-direction, also running variable for counting average data output
@@ -183,10 +234,20 @@
       INTEGER(iwp) ::  l                 !< running index for surface-element orientation
       INTEGER(iwp) ::  m                 !< running index for surface elements
+      INTEGER(iwp) ::  mm                !< local counting variable for surface elements
       INTEGER(iwp) ::  npg               !< counter variable for all surface elements ( or polygons )
       INTEGER(iwp) ::  point_index_count !< local counter variable for point index
+      INTEGER(iwp) ::  start_count       !< local start counter for the surface index
       
       INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_points_on_pe   !< array which contains the number of points on all mpi ranks
+      INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_surfaces_on_pe !< array which contains the number of surfaces on all mpi ranks
       INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) ::  point_index !< dummy array used to check where the reference points for surface polygons are located
            
+      REAL(wp) ::  az    !< azimuth angle, indicated the vertical orientation of a surface element
+      REAL(wp) ::  off_x !< grid offset in x-direction between the stored grid index and the actual wall
+      REAL(wp) ::  off_y !< grid offset in y-direction between the stored grid index and the actual wall
+
+      REAL(wp), DIMENSION(:), ALLOCATABLE ::  netcdf_data_1d  !< dummy array to output 1D data into netcdf file
+
+
 !
 !--   Determine the number of surface elements on subdomain
@@ -198,4 +259,12 @@
                   + surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns       !eastward-facing   
 !
+!--    Determine the total number of surfaces in the model domain
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1,                  &
+                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
+#else
+       surfaces%ns_total = surfaces%ns
+#endif
+!
 !--   Allocate output variable and set to _FillValue attribute
       ALLOCATE ( surfaces%var_out(1:surfaces%ns) )
@@ -209,22 +278,46 @@
       ENDIF
 !
-!--   Initialize points and polygon data. 
+!--   If output to VTK format is enabled, initialize point and polygon data. 
 !--   In a first step, count the number of points which are defining 
 !--   the surfaces and the polygons.
-      ALLOCATE( point_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
-      point_index = -1
-!
-!--   Horizontal default surfaces
-      surfaces%npoints = 0
-      DO  l = 0, 1
-         DO  m = 1, surf_def_h(0)%ns
-!
-!--         Determine the indices of the respective grid cell inside the topography
-            i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
-            j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
-            k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
-!
-!--         Check if the vertices that define the surface element are already 
-!--         defined, if not, increment the counter.
+      IF ( to_vtk )  THEN
+         ALLOCATE( point_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+         point_index = -1
+!
+!--      Horizontal default surfaces
+         surfaces%npoints = 0
+         DO  l = 0, 1
+            DO  m = 1, surf_def_h(0)%ns
+!
+!--            Determine the indices of the respective grid cell inside the topography
+               i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
+               j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
+               k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
+!
+!--            Check if the vertices that define the surface element are already 
+!--            defined, if not, increment the counter.
+               IF ( point_index(k,j,i) < 0 )  THEN
+                  surfaces%npoints   = surfaces%npoints + 1    
+                  point_index(k,j,i) = surfaces%npoints - 1
+               ENDIF
+               IF ( point_index(k,j,i+1) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j,i+1) = surfaces%npoints - 1
+               ENDIF
+               IF ( point_index(k,j+1,i+1) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k,j+1,i+1) = surfaces%npoints - 1
+               ENDIF
+               IF ( point_index(k,j+1,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j+1,i) = surfaces%npoints - 1
+               ENDIF           
+            ENDDO
+         ENDDO
+         DO  m = 1, surf_lsm_h%ns
+            i = surf_lsm_h%i(m) + surf_lsm_h%ioff
+            j = surf_lsm_h%j(m) + surf_lsm_h%joff
+            k = surf_lsm_h%k(m) + surf_lsm_h%koff
+
             IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
@@ -242,231 +335,253 @@
                surfaces%npoints     = surfaces%npoints + 1    
                point_index(k,j+1,i) = surfaces%npoints - 1
-            ENDIF           
+            ENDIF    
          ENDDO
-      ENDDO
-      DO  m = 1, surf_lsm_h%ns
-         i = surf_lsm_h%i(m) + surf_lsm_h%ioff
-         j = surf_lsm_h%j(m) + surf_lsm_h%joff
-         k = surf_lsm_h%k(m) + surf_lsm_h%koff
-
-         IF ( point_index(k,j,i) < 0 )  THEN
-            surfaces%npoints   = surfaces%npoints + 1    
-            point_index(k,j,i) = surfaces%npoints - 1
-         ENDIF
-         IF ( point_index(k,j,i+1) < 0 )  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j,i+1) = surfaces%npoints - 1
-         ENDIF
-         IF ( point_index(k,j+1,i+1) < 0 )  THEN
-            surfaces%npoints       = surfaces%npoints + 1    
-            point_index(k,j+1,i+1) = surfaces%npoints - 1
-         ENDIF
-         IF ( point_index(k,j+1,i) < 0 )  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j+1,i) = surfaces%npoints - 1
-         ENDIF    
-      ENDDO
-      DO  m = 1, surf_usm_h%ns
-         i = surf_usm_h%i(m) + surf_usm_h%ioff
-         j = surf_usm_h%j(m) + surf_usm_h%joff
-         k = surf_usm_h%k(m) + surf_usm_h%koff
-
-         IF ( point_index(k,j,i) < 0 )  THEN
-            surfaces%npoints   = surfaces%npoints + 1    
-            point_index(k,j,i) = surfaces%npoints - 1
-         ENDIF
-         IF ( point_index(k,j,i+1) < 0 )  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j,i+1) = surfaces%npoints - 1
-         ENDIF
-         IF ( point_index(k,j+1,i+1) < 0 )  THEN
-            surfaces%npoints       = surfaces%npoints + 1    
-            point_index(k,j+1,i+1) = surfaces%npoints - 1
-         ENDIF
-         IF ( point_index(k,j+1,i) < 0 )  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j+1,i) = surfaces%npoints - 1
-         ENDIF     
-      ENDDO
-!
-!--   Vertical surfaces
-      DO  l = 0, 3
-         DO  m = 1, surf_def_v(l)%ns
-!
-!--         Determine the indices of the respective grid cell inside the topography.
-!--         Please note, j-index for north-facing surfaces ( l==0 ) is 
-!--         identical to the reference j-index outside the grid box. 
-!--         Equivalent for east-facing surfaces and i-index.
-            i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
-            j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
-            k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
-!
-!--         Lower left /front vertex
+         DO  m = 1, surf_usm_h%ns
+            i = surf_usm_h%i(m) + surf_usm_h%ioff
+            j = surf_usm_h%j(m) + surf_usm_h%joff
+            k = surf_usm_h%k(m) + surf_usm_h%koff
+
             IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
                point_index(k,j,i) = surfaces%npoints - 1
-            ENDIF 
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .OR.  l == 1 )  THEN
+            ENDIF
+            IF ( point_index(k,j,i+1) < 0 )  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k,j,i+1) = surfaces%npoints - 1
+            ENDIF
+            IF ( point_index(k,j+1,i+1) < 0 )  THEN
+               surfaces%npoints       = surfaces%npoints + 1    
+               point_index(k,j+1,i+1) = surfaces%npoints - 1
+            ENDIF
+            IF ( point_index(k,j+1,i) < 0 )  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k,j+1,i) = surfaces%npoints - 1
+            ENDIF     
+         ENDDO
+!
+!--      Vertical surfaces
+         DO  l = 0, 3
+            DO  m = 1, surf_def_v(l)%ns
+!
+!--            Determine the indices of the respective grid cell inside the
+!--            topography. Please note, j-index for north-facing surfaces 
+!--            ( l==0 ) is identical to the reference j-index outside the grid  
+!--            box. Equivalent for east-facing surfaces and i-index.
+               i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
+               j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
+               k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
+!
+!--            Lower left /front vertex
+               IF ( point_index(k,j,i) < 0 )  THEN
+                  surfaces%npoints   = surfaces%npoints + 1    
+                  point_index(k,j,i) = surfaces%npoints - 1
+               ENDIF 
+!
+!--            Upper / lower right index for north- and south-facing surfaces
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  IF ( point_index(k,j,i+1) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j,i+1) = surfaces%npoints - 1
+                  ENDIF 
+                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j,i+1) = surfaces%npoints - 1
+                  ENDIF
+!
+!--            Upper / lower front index for east- and west-facing surfaces
+               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+                  IF ( point_index(k,j+1,i) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j+1,i) = surfaces%npoints - 1
+                  ENDIF 
+                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j+1,i) = surfaces%npoints - 1
+                  ENDIF
+               ENDIF
+!
+!--            Upper left / front vertex
+               IF ( point_index(k+1,j,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k+1,j,i) = surfaces%npoints - 1
+               ENDIF
+            ENDDO
+            DO  m = 1, surf_lsm_v(l)%ns
+               i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
+               j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
+               k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
+!
+!--            Lower left /front vertex
+               IF ( point_index(k,j,i) < 0 )  THEN
+                  surfaces%npoints   = surfaces%npoints + 1    
+                  point_index(k,j,i) = surfaces%npoints - 1
+               ENDIF 
+!
+!--            Upper / lower right index for north- and south-facing surfaces
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  IF ( point_index(k,j,i+1) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j,i+1) = surfaces%npoints - 1
+                  ENDIF 
+                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j,i+1) = surfaces%npoints - 1
+                  ENDIF
+!
+!--            Upper / lower front index for east- and west-facing surfaces
+               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+                  IF ( point_index(k,j+1,i) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j+1,i) = surfaces%npoints - 1
+                  ENDIF 
+                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j+1,i) = surfaces%npoints - 1
+                  ENDIF
+               ENDIF
+!
+!--            Upper left / front vertex
+               IF ( point_index(k+1,j,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k+1,j,i) = surfaces%npoints - 1
+               ENDIF
+            ENDDO
+
+            DO  m = 1, surf_usm_v(l)%ns
+               i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
+               j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
+               k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
+!
+!--            Lower left /front vertex
+               IF ( point_index(k,j,i) < 0 )  THEN
+                  surfaces%npoints   = surfaces%npoints + 1    
+                  point_index(k,j,i) = surfaces%npoints - 1
+               ENDIF 
+!
+!--            Upper / lower right index for north- and south-facing surfaces
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  IF ( point_index(k,j,i+1) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j,i+1) = surfaces%npoints - 1
+                  ENDIF 
+                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j,i+1) = surfaces%npoints - 1
+                  ENDIF
+!
+!--            Upper / lower front index for east- and west-facing surfaces
+               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+                  IF ( point_index(k,j+1,i) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j+1,i) = surfaces%npoints - 1
+                  ENDIF 
+                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j+1,i) = surfaces%npoints - 1
+                  ENDIF
+               ENDIF
+!
+!--            Upper left / front vertex
+               IF ( point_index(k+1,j,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k+1,j,i) = surfaces%npoints - 1
+               ENDIF
+            ENDDO
+
+         ENDDO
+!
+!--      Allocate the number of points and polygons. Note, the number of polygons
+!--      is identical to the number of surfaces elements, whereas the number
+!--      of points (vertices), which define the polygons, can be larger.
+         ALLOCATE( surfaces%points(3,1:surfaces%npoints) )
+         ALLOCATE( surfaces%polygons(5,1:surfaces%ns)    )
+!
+!--      Note, PARAVIEW expects consecutively ordered points, in order to
+!--      unambiguously identify surfaces. 
+!--      Hence, all PEs should know where they start counting, depending on the
+!--      number of points on the other PE's with lower MPI rank. 
+#if defined( __parallel )
+         CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER,                 &
+                             num_points_on_pe, 1, MPI_INTEGER, comm2d, ierr  )
+#else
+         num_points_on_pe = surfaces%npoints
+#endif
+
+!
+!--      After the number of vertices is counted, repeat the loops and define the
+!--      vertices. Start with the horizontal default surfaces.
+!--      First, however, determine the offset where couting of points should be
+!--      started, which is the sum of points of all PE's with lower MPI rank.
+         i                 = 0
+         point_index_count = 0
+         DO WHILE ( i < myid  .AND.  i <= SIZE( num_points_on_pe ) )
+            point_index_count = point_index_count + num_points_on_pe(i)
+            i                 = i + 1
+         ENDDO
+          
+         surfaces%npoints = 0
+         point_index      = -1
+         npg              = 0
+      
+         DO  l = 0, 1
+            DO  m = 1, surf_def_h(0)%ns
+!
+!--            Determine the indices of the respective grid cell inside the 
+!--            topography.
+               i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
+               j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
+               k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
+!
+!--            Check if the vertices that define the surface element are  
+!--            already defined, if not, increment the counter.
+               IF ( point_index(k,j,i) < 0 )  THEN
+                  surfaces%npoints   = surfaces%npoints + 1    
+                  point_index(k,j,i) = point_index_count
+                  point_index_count  = point_index_count + 1               
+                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
                IF ( point_index(k,j,i+1) < 0 )  THEN
                   surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints - 1
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                  point_index(k,j,i+1) = point_index_count
+                  point_index_count    = point_index_count + 1 
+                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
+               IF ( point_index(k,j+1,i+1) < 0 )  THEN
                   surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints - 1
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+                  point_index(k,j+1,i+1) = point_index_count
+                  point_index_count      = point_index_count + 1 
+                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
                IF ( point_index(k,j+1,i) < 0 )  THEN
                   surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints - 1
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints - 1
-               ENDIF
-            ENDIF
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints - 1
-            ENDIF
+                  point_index(k,j+1,i) = point_index_count
+                  point_index_count    = point_index_count + 1 
+                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
+               
+               npg                        = npg + 1
+               surfaces%polygons(1,npg)   = 4
+               surfaces%polygons(2,npg)   = point_index(k,j,i)
+               surfaces%polygons(3,npg)   = point_index(k,j,i+1)
+               surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
+               surfaces%polygons(5,npg)   = point_index(k,j+1,i)
+            ENDDO
          ENDDO
-         DO  m = 1, surf_lsm_v(l)%ns
-            i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
-            j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
-            k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
-!
-!--         Lower left /front vertex
-            IF ( point_index(k,j,i) < 0 )  THEN
-               surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints - 1
-            ENDIF 
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .OR.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints - 1
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints - 1
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints - 1
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints - 1
-               ENDIF
-            ENDIF
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints - 1
-            ENDIF
-         ENDDO
-
-         DO  m = 1, surf_usm_v(l)%ns
-            i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
-            j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
-            k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
-!
-!--         Lower left /front vertex
-            IF ( point_index(k,j,i) < 0 )  THEN
-               surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints - 1
-            ENDIF 
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .OR.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints - 1
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints - 1
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints - 1
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints - 1
-               ENDIF
-            ENDIF
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints - 1
-            ENDIF
-         ENDDO
-
-      ENDDO
-!
-!--   Allocate the number of points and polygons. Note, the number of polygons
-!--   is identical to the number of surfaces elements, whereas the number
-!--   of points (vertices) which define the polygons can be larger.
-      ALLOCATE( surfaces%points(3,1:surfaces%npoints) )
-      ALLOCATE( surfaces%polygons(5,1:surfaces%ns)    )
-!
-!--   Note, PARAVIEW expects consecutively ordered points which can be
-!--   unambiguously identified. 
-!--   Hence, all PEs know where they should start counting, depending on the
-!--   number of points on the other PE's with lower MPI rank. 
-#if defined( __parallel )
-      CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER,                    &
-                          num_points_on_pe, 1, MPI_INTEGER, comm2d, ierr  )
-#else
-      num_points_on_pe = surfaces%npoints
-#endif
-
-!
-!--   After the number of vertices is counted, repeat the loops and define the
-!--   vertices. Start with the horizontal default surfaces
-!--   First, however, determine the offset where couting of points should be
-!--   started, which is the sum of points of all PE's with lower MPI rank.
-      i                 = 0
-      point_index_count = 0
-      DO WHILE ( i < myid  .AND.  i <= SIZE( num_points_on_pe ) )
-         point_index_count = point_index_count + num_points_on_pe(i)
-         i                 = i + 1
-      ENDDO
-          
-      
-      surfaces%npoints = 0
-      point_index      = -1
-      npg              = 0
-      
-      DO  l = 0, 1
-         DO  m = 1, surf_def_h(0)%ns
-!
-!--         Determine the indices of the respective grid cell inside the topography
-            i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
-            j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
-            k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
-!
-!--         Check if the vertices that define the surface element are already 
-!--         defined, if not, increment the counter.
+         DO  m = 1, surf_lsm_h%ns
+            i = surf_lsm_h%i(m) + surf_lsm_h%ioff
+            j = surf_lsm_h%j(m) + surf_lsm_h%joff
+            k = surf_lsm_h%k(m) + surf_lsm_h%koff
             IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
                point_index(k,j,i) = point_index_count
-               point_index_count  = point_index_count + 1               
+               point_index_count  = point_index_count + 1 
                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
@@ -503,110 +618,12 @@
             surfaces%polygons(3,npg)   = point_index(k,j,i+1)
             surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
-            surfaces%polygons(5,npg)   = point_index(k,j+1,i)
+            surfaces%polygons(5,npg)   = point_index(k,j+1,i) 
          ENDDO
-      ENDDO
-      DO  m = 1, surf_lsm_h%ns
-         i = surf_lsm_h%i(m) + surf_lsm_h%ioff
-         j = surf_lsm_h%j(m) + surf_lsm_h%joff
-         k = surf_lsm_h%k(m) + surf_lsm_h%koff
-         IF ( point_index(k,j,i) < 0 )  THEN
-            surfaces%npoints   = surfaces%npoints + 1    
-            point_index(k,j,i) = point_index_count
-            point_index_count  = point_index_count + 1 
-            surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         IF ( point_index(k,j,i+1) < 0 )  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j,i+1) = point_index_count
-            point_index_count    = point_index_count + 1 
-            surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         IF ( point_index(k,j+1,i+1) < 0 )  THEN
-            surfaces%npoints       = surfaces%npoints + 1    
-            point_index(k,j+1,i+1) = point_index_count
-            point_index_count      = point_index_count + 1 
-            surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         IF ( point_index(k,j+1,i) < 0 )  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j+1,i) = point_index_count
-            point_index_count    = point_index_count + 1 
-            surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         
-         npg                        = npg + 1
-         surfaces%polygons(1,npg)   = 4
-         surfaces%polygons(2,npg)   = point_index(k,j,i)
-         surfaces%polygons(3,npg)   = point_index(k,j,i+1)
-         surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
-         surfaces%polygons(5,npg)   = point_index(k,j+1,i) 
-      ENDDO
-
-      DO  m = 1, surf_usm_h%ns
-         i = surf_usm_h%i(m) + surf_usm_h%ioff
-         j = surf_usm_h%j(m) + surf_usm_h%joff
-         k = surf_usm_h%k(m) + surf_usm_h%koff
-
-         IF ( point_index(k,j,i) < 0 )  THEN
-            surfaces%npoints   = surfaces%npoints + 1    
-            point_index(k,j,i) = point_index_count
-            point_index_count  = point_index_count + 1 
-            surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         IF ( point_index(k,j,i+1) < 0 )  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j,i+1) = point_index_count
-            point_index_count    = point_index_count + 1 
-            surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         IF ( point_index(k,j+1,i+1) < 0 )  THEN
-            surfaces%npoints       = surfaces%npoints + 1    
-            point_index(k,j+1,i+1) = point_index_count
-            point_index_count      = point_index_count + 1 
-            surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         IF ( point_index(k,j+1,i) < 0 )  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j+1,i) = point_index_count
-            point_index_count    = point_index_count + 1 
-            surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         
-         npg                        = npg + 1
-         surfaces%polygons(1,npg)   = 4
-         surfaces%polygons(2,npg)   = point_index(k,j,i)
-         surfaces%polygons(3,npg)   = point_index(k,j,i+1)
-         surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
-         surfaces%polygons(5,npg)   = point_index(k,j+1,i)    
-      ENDDO
-
-      DO  l = 0, 3
-         DO  m = 1, surf_def_v(l)%ns
-!
-!--         Determine the indices of the respective grid cell inside the topography.
-!--         Please note, j-index for north-facing surfaces ( l==0 ) is 
-!--         identical to the reference j-index outside the grid box. 
-!--         Equivalent for east-facing surfaces and i-index.
-            i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
-            j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
-            k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
-!
-!--         Lower left /front vertex
+     
+         DO  m = 1, surf_usm_h%ns
+            i = surf_usm_h%i(m) + surf_usm_h%ioff
+            j = surf_usm_h%j(m) + surf_usm_h%joff
+            k = surf_usm_h%k(m) + surf_usm_h%koff
+     
             IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
@@ -615,249 +632,639 @@
                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k-1)
-            ENDIF 
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .OR.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = point_index_count
-                  point_index_count    = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k-1)
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = point_index_count
-                  point_index_count      = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = point_index_count
-                  point_index_count    = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k-1)
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = point_index_count
-                  point_index_count      = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF
+               surfaces%points(3,surfaces%npoints) = zw(k)
             ENDIF
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0 )  THEN
+            IF ( point_index(k,j,i+1) < 0 )  THEN
                surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = point_index_count
-               point_index_count    = point_index_count + 1
-               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
+               point_index(k,j,i+1) = point_index_count
+               point_index_count    = point_index_count + 1 
+               surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+               surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+               surfaces%points(3,surfaces%npoints) = zw(k)
+            ENDIF
+            IF ( point_index(k,j+1,i+1) < 0 )  THEN
+               surfaces%npoints       = surfaces%npoints + 1    
+               point_index(k,j+1,i+1) = point_index_count
+               point_index_count      = point_index_count + 1 
+               surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+               surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+               surfaces%points(3,surfaces%npoints) = zw(k)
+            ENDIF
+            IF ( point_index(k,j+1,i) < 0 )  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k,j+1,i) = point_index_count
+               point_index_count    = point_index_count + 1 
+               surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+               surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
                surfaces%points(3,surfaces%npoints) = zw(k)
             ENDIF
             
-            npg = npg + 1            
-            IF ( l == 0  .OR.  l == 1 )  THEN
-               surfaces%polygons(1,npg)   = 4
-               surfaces%polygons(2,npg)   = point_index(k,j,i)
-               surfaces%polygons(3,npg)   = point_index(k,j,i+1)
-               surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
-               surfaces%polygons(5,npg)   = point_index(k+1,j,i)  
-            ELSE
-               surfaces%polygons(1,npg)   = 4
-               surfaces%polygons(2,npg)   = point_index(k,j,i)
-               surfaces%polygons(3,npg)   = point_index(k,j+1,i)
-               surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
-               surfaces%polygons(5,npg)   = point_index(k+1,j,i)
-            ENDIF
-            
+            npg                        = npg + 1
+            surfaces%polygons(1,npg)   = 4
+            surfaces%polygons(2,npg)   = point_index(k,j,i)
+            surfaces%polygons(3,npg)   = point_index(k,j,i+1)
+            surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
+            surfaces%polygons(5,npg)   = point_index(k,j+1,i)    
          ENDDO
 
-         DO  m = 1, surf_lsm_v(l)%ns
-            i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
-            j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
-            k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
-!
-!--         Lower left /front vertex
-            IF ( point_index(k,j,i) < 0 )  THEN
-               surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = point_index_count
-               point_index_count  = point_index_count + 1
-               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k-1)
-            ENDIF 
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .OR.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = point_index_count
-                  point_index_count    = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+         DO  l = 0, 3
+            DO  m = 1, surf_def_v(l)%ns
+!        
+!--            Determine the indices of the respective grid cell inside the topography.
+!--            Please note, j-index for north-facing surfaces ( l==0 ) is 
+!--            identical to the reference j-index outside the grid box. 
+!--            Equivalent for east-facing surfaces and i-index.
+               i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
+               j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
+               k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
+!        
+!--            Lower left /front vertex
+               IF ( point_index(k,j,i) < 0 )  THEN
+                  surfaces%npoints   = surfaces%npoints + 1    
+                  point_index(k,j,i) = point_index_count
+                  point_index_count  = point_index_count + 1 
+                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k-1)
                ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = point_index_count
-                  point_index_count      = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+!        
+!--            Upper / lower right index for north- and south-facing surfaces
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  IF ( point_index(k,j,i+1) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j,i+1) = point_index_count
+                     point_index_count    = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k-1)
+                  ENDIF 
+                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j,i+1) = point_index_count
+                     point_index_count      = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k)
+                  ENDIF
+!        
+!--            Upper / lower front index for east- and west-facing surfaces
+               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+                  IF ( point_index(k,j+1,i) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j+1,i) = point_index_count
+                     point_index_count    = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k-1)
+                  ENDIF 
+                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j+1,i) = point_index_count
+                     point_index_count      = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k)
+                  ENDIF
+               ENDIF
+!        
+!--            Upper left / front vertex
+               IF ( point_index(k+1,j,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k+1,j,i) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k)
                ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = point_index_count
-                  point_index_count    = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+               
+               npg = npg + 1            
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  surfaces%polygons(1,npg)   = 4
+                  surfaces%polygons(2,npg)   = point_index(k,j,i)
+                  surfaces%polygons(3,npg)   = point_index(k,j,i+1)
+                  surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)  
+               ELSE
+                  surfaces%polygons(1,npg)   = 4
+                  surfaces%polygons(2,npg)   = point_index(k,j,i)
+                  surfaces%polygons(3,npg)   = point_index(k,j+1,i)
+                  surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
+               ENDIF
+               
+            ENDDO
+         
+            DO  m = 1, surf_lsm_v(l)%ns
+               i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
+               j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
+               k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
+!        
+!--            Lower left /front vertex
+               IF ( point_index(k,j,i) < 0 )  THEN
+                  surfaces%npoints   = surfaces%npoints + 1    
+                  point_index(k,j,i) = point_index_count
+                  point_index_count  = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k-1)
                ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = point_index_count
-                  point_index_count      = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+!        
+!--            Upper / lower right index for north- and south-facing surfaces
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  IF ( point_index(k,j,i+1) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j,i+1) = point_index_count
+                     point_index_count    = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k-1)
+                  ENDIF 
+                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j,i+1) = point_index_count
+                     point_index_count      = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k)
+                  ENDIF
+!        
+!--            Upper / lower front index for east- and west-facing surfaces
+               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+                  IF ( point_index(k,j+1,i) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j+1,i) = point_index_count
+                     point_index_count    = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k-1)
+                  ENDIF 
+                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j+1,i) = point_index_count
+                     point_index_count      = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k)
+                  ENDIF
+               ENDIF
+!        
+!--            Upper left / front vertex
+               IF ( point_index(k+1,j,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k+1,j,i) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k)
                ENDIF
-            ENDIF
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = point_index_count
-               point_index_count    = point_index_count + 1
-               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k)
-            ENDIF
-            
-            npg = npg + 1            
-            IF ( l == 0  .OR.  l == 1 )  THEN
-               surfaces%polygons(1,npg)   = 4
-               surfaces%polygons(2,npg)   = point_index(k,j,i)
-               surfaces%polygons(3,npg)   = point_index(k,j,i+1)
-               surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
-               surfaces%polygons(5,npg)   = point_index(k+1,j,i)  
-            ELSE
-               surfaces%polygons(1,npg)   = 4
-               surfaces%polygons(2,npg)   = point_index(k,j,i)
-               surfaces%polygons(3,npg)   = point_index(k,j+1,i)
-               surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
-               surfaces%polygons(5,npg)   = point_index(k+1,j,i)
-            ENDIF
-         ENDDO
-         DO  m = 1, surf_usm_v(l)%ns
-            i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
-            j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
-            k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
-!
-!--         Lower left /front vertex
-            IF ( point_index(k,j,i) < 0 )  THEN
-               surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = point_index_count
-               point_index_count  = point_index_count + 1
-               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k-1)
-            ENDIF 
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .OR.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = point_index_count
-                  point_index_count    = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+               
+               npg = npg + 1            
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  surfaces%polygons(1,npg)   = 4
+                  surfaces%polygons(2,npg)   = point_index(k,j,i)
+                  surfaces%polygons(3,npg)   = point_index(k,j,i+1)
+                  surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)  
+               ELSE
+                  surfaces%polygons(1,npg)   = 4
+                  surfaces%polygons(2,npg)   = point_index(k,j,i)
+                  surfaces%polygons(3,npg)   = point_index(k,j+1,i)
+                  surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
+               ENDIF
+            ENDDO
+            DO  m = 1, surf_usm_v(l)%ns
+               i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
+               j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
+               k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
+!        
+!--            Lower left /front vertex
+               IF ( point_index(k,j,i) < 0 )  THEN
+                  surfaces%npoints   = surfaces%npoints + 1    
+                  point_index(k,j,i) = point_index_count
+                  point_index_count  = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k-1)
                ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = point_index_count
-                  point_index_count      = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+!        
+!--            Upper / lower right index for north- and south-facing surfaces
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  IF ( point_index(k,j,i+1) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j,i+1) = point_index_count
+                     point_index_count    = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k-1)
+                  ENDIF 
+                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j,i+1) = point_index_count
+                     point_index_count      = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k)
+                  ENDIF
+!        
+!--            Upper / lower front index for east- and west-facing surfaces
+               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+                  IF ( point_index(k,j+1,i) < 0 )  THEN
+                     surfaces%npoints     = surfaces%npoints + 1    
+                     point_index(k,j+1,i) = point_index_count
+                     point_index_count    = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k-1)
+                  ENDIF 
+                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                     surfaces%npoints       = surfaces%npoints + 1    
+                     point_index(k+1,j+1,i) = point_index_count
+                     point_index_count      = point_index_count + 1
+                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                     surfaces%points(3,surfaces%npoints) = zw(k)
+                  ENDIF
+               ENDIF
+!        
+!--            Upper left / front vertex
+               IF ( point_index(k+1,j,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k+1,j,i) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k)
                ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = point_index_count
-                  point_index_count    = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k-1)
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = point_index_count
-                  point_index_count      = point_index_count + 1
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF
-            ENDIF
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = point_index_count
-               point_index_count    = point_index_count + 1
-               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k)
-            ENDIF
-            
-            npg = npg + 1            
-            IF ( l == 0  .OR.  l == 1 )  THEN
-               surfaces%polygons(1,npg)   = 4
-               surfaces%polygons(2,npg)   = point_index(k,j,i)
-               surfaces%polygons(3,npg)   = point_index(k,j,i+1)
-               surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
-               surfaces%polygons(5,npg)   = point_index(k+1,j,i)  
-            ELSE
-               surfaces%polygons(1,npg)   = 4
-               surfaces%polygons(2,npg)   = point_index(k,j,i)
-               surfaces%polygons(3,npg)   = point_index(k,j+1,i)
-               surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
-               surfaces%polygons(5,npg)   = point_index(k+1,j,i)
-            ENDIF
+               
+               npg = npg + 1            
+               IF ( l == 0  .OR.  l == 1 )  THEN
+                  surfaces%polygons(1,npg)   = 4
+                  surfaces%polygons(2,npg)   = point_index(k,j,i)
+                  surfaces%polygons(3,npg)   = point_index(k,j,i+1)
+                  surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)  
+               ELSE
+                  surfaces%polygons(1,npg)   = 4
+                  surfaces%polygons(2,npg)   = point_index(k,j,i)
+                  surfaces%polygons(3,npg)   = point_index(k,j+1,i)
+                  surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
+               ENDIF
+            ENDDO
+         
          ENDDO
-
-      ENDDO
-!
-!--   Deallocate temporary dummy variable
-      DEALLOCATE ( point_index )
-!
-!--   Sum-up total number of surface elements and vertices on domain. This
-!--   will be needed for post-processing. 
-      surfaces%npoints_total = 0
+!        
+!--      Deallocate temporary dummy variable
+         DEALLOCATE ( point_index )
+!
+!--      Sum-up total number of vertices on domain. This
+!--      will be needed for post-processing. 
+         surfaces%npoints_total = 0
 #if defined( __parallel )
-       CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1,        &
-                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
-       CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1,                  &
-                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
+          CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1,     &
+                              MPI_INTEGER, MPI_SUM, comm2d, ierr )
 #else
-       surfaces%npoints_total = surfaces%npoints
-       surfaces%ns_total      = surfaces%ns
+          surfaces%npoints_total = surfaces%npoints
 #endif
+       ENDIF
+!
+!--    If output to netcdf is enabled, set-up the coordinate arrays that 
+!--    unambiguously describe the position and orientation of each surface 
+!--    element. 
+       IF ( to_netcdf )  THEN
+!
+!--       Allocate local coordinate arrays
+          ALLOCATE( surfaces%s(1:surfaces%ns)       )
+          ALLOCATE( surfaces%xs(1:surfaces%ns)      )
+          ALLOCATE( surfaces%ys(1:surfaces%ns)      )
+          ALLOCATE( surfaces%zs(1:surfaces%ns)      )
+          ALLOCATE( surfaces%azimuth(1:surfaces%ns) )
+          ALLOCATE( surfaces%zenith(1:surfaces%ns)  )
+          ALLOCATE( surfaces%es_utm(1:surfaces%ns)  )
+          ALLOCATE( surfaces%ns_utm(1:surfaces%ns)  )
+!
+!--       Gather the number of surface on each processor, in order to number
+!--       the surface elements in ascending order with respect to the total 
+!--       number of surfaces in the domain.
+#if defined( __parallel )
+          CALL MPI_ALLGATHER( surfaces%ns, 1, MPI_INTEGER,                     &
+                              num_surfaces_on_pe, 1, MPI_INTEGER, comm2d, ierr  )
+#else
+          num_surfaces_on_pe = surfaces%ns
+#endif
+!
+!--       First, however, determine the offset where couting of the surfaces 
+!--       should start (the sum of surfaces on all PE's with lower MPI rank).
+          i           = 0
+          start_count = 1
+          DO WHILE ( i < myid  .AND.  i <= SIZE( num_surfaces_on_pe ) )
+             start_count = start_count + num_surfaces_on_pe(i)
+             i           = i + 1
+          ENDDO
+!       
+!--       Set coordinate arrays. For horizontal surfaces, azimuth 
+!--       angles are not defined (fill value). Zenith angle is 0 (180) for 
+!--       upward (downward)-facing surfaces.
+          i  = start_count
+          mm = 1
+          DO  m = 1, surf_def_h(0)%ns
+             surfaces%s(mm)       = i
+             surfaces%xs(mm)      = ( surf_def_h(0)%i(m) + 0.5_wp ) * dx
+             surfaces%ys(mm)      = ( surf_def_h(0)%j(m) + 0.5_wp ) * dy
+             surfaces%zs(mm)      = zw(surf_def_h(0)%k(m)+surf_def_h(0)%koff)
+             surfaces%azimuth(mm) = surfaces%fillvalue
+             surfaces%zenith(mm)  = 0.0
+             i                    = i + 1
+             mm                   = mm + 1
+          ENDDO
+          DO  m = 1, surf_lsm_h%ns
+             surfaces%s(mm)       = i
+             surfaces%xs(mm)      = ( surf_lsm_h%i(m) + 0.5_wp ) * dx
+             surfaces%ys(mm)      = ( surf_lsm_h%j(m) + 0.5_wp ) * dy
+             surfaces%zs(mm)      = zw(surf_lsm_h%k(m)+surf_lsm_h%koff)
+             surfaces%azimuth(mm) = surfaces%fillvalue
+             surfaces%zenith(mm)  = 0.0
+             i                    = i + 1
+             mm                   = mm + 1
+          ENDDO
+          DO  m = 1, surf_usm_h%ns
+             surfaces%s(mm)       = i
+             surfaces%xs(mm)      = ( surf_usm_h%i(m) + 0.5_wp ) * dx
+             surfaces%ys(mm)      = ( surf_usm_h%j(m) + 0.5_wp ) * dy
+             surfaces%zs(mm)      = zw(surf_usm_h%k(m)+surf_usm_h%koff)
+             surfaces%azimuth(mm) = surfaces%fillvalue
+             surfaces%zenith(mm)  = 0.0
+             i                    = i + 1
+             mm                   = mm + 1
+          ENDDO
+          DO  m = 1, surf_def_h(1)%ns
+             surfaces%s(mm)       = i
+             surfaces%xs(mm)      = ( surf_def_h(1)%i(m) + 0.5_wp ) * dx
+             surfaces%ys(mm)      = ( surf_def_h(1)%j(m) + 0.5_wp ) * dy
+             surfaces%zs(mm)      = zw(surf_def_h(1)%k(m)+surf_def_h(1)%koff)
+             surfaces%azimuth(mm) = surfaces%fillvalue
+             surfaces%zenith(mm)  = 180.0
+             i                    = i + 1
+             mm                   = mm + 1
+          ENDDO
+!       
+!--       For vertical surfaces, zenith angles are not defined (fill value).  
+!--       Azimuth angle: eastward (0), westward (180), northward (270), 
+!--       southward (90). 
+          DO  l = 0, 3
+             IF ( l == 0 )  THEN
+                az    = 270.0_wp
+                off_x =   0.0_wp
+                off_y = - 0.5_wp
+             ELSEIF ( l == 1 )  THEN
+                az    = 90.0_wp
+                off_x =  0.0_wp
+                off_y =  0.5_wp
+             ELSEIF ( l == 2 )  THEN
+                az    =   0.0_wp
+                off_x = - 0.5_wp
+                off_y =   0.0_wp
+             ELSEIF ( l == 3 )  THEN
+                az   = 180.0_wp
+                off_x =  0.5_wp
+                off_y =  0.0_wp
+             ENDIF
+                
+             DO  m = 1, surf_def_v(l)%ns
+                surfaces%s(mm)       = i
+                surfaces%xs(mm)      = ( surf_def_v(l)%i(m) + off_x ) * dx
+                surfaces%ys(mm)      = ( surf_def_v(l)%j(m) + off_y ) * dy
+                surfaces%zs(mm)      = zu(surf_def_v(l)%k(m))
+                surfaces%azimuth(mm) = az
+                surfaces%zenith(mm)  = surfaces%fillvalue
+                i                    = i + 1
+                mm                   = mm + 1
+             ENDDO
+             DO  m = 1, surf_lsm_v(l)%ns
+                surfaces%s(mm)       = i
+                surfaces%xs(mm)      = ( surf_lsm_v(l)%i(m) + off_x ) * dx
+                surfaces%ys(mm)      = ( surf_lsm_v(l)%j(m) + off_y ) * dy
+                surfaces%zs(mm)      = zu(surf_lsm_v(l)%k(m))
+                surfaces%azimuth(mm) = az
+                surfaces%zenith(mm)  = surfaces%fillvalue
+                i                    = i + 1
+                mm                   = mm + 1
+             ENDDO
+             DO  m = 1, surf_usm_v(l)%ns
+                surfaces%s(mm)       = i
+                surfaces%xs(mm)      = ( surf_usm_v(l)%i(m) + off_x ) * dx
+                surfaces%ys(mm)      = ( surf_usm_v(l)%j(m) + off_y ) * dy
+                surfaces%zs(mm)      = zu(surf_usm_v(l)%k(m))
+                surfaces%azimuth(mm) = az
+                surfaces%zenith(mm)  = surfaces%fillvalue
+                i                    = i + 1
+                mm                   = mm + 1
+             ENDDO
+          ENDDO
+!       
+!--       Finally, define UTM coordinates, which are the x/y-coordinates 
+!--       plus the origin (lower-left coordinate of the model domain).
+          surfaces%es_utm = surfaces%xs + init_model%origin_x 
+          surfaces%ns_utm = surfaces%ys + init_model%origin_y 
+!
+!--       Initialize NetCDF data output. Please note, local start position for 
+!--       the surface elements in the NetCDF file is surfaces%s(1), while 
+!--       the number of surfaces on the subdomain is given by surfaces%ns.
+#if defined( __netcdf4_parallel )
+
+!
+!--       Calculate number of time steps to be output
+          ntdim_surf(0) = dosurf_time_count(0) + CEILING(                            &
+                        ( end_time - MAX(                                            &
+                            MERGE(skip_time_dosurf, skip_time_dosurf + spinup_time,  &
+                                  data_output_during_spinup ),                       &
+                            simulated_time_at_begin )                                &
+                        ) / dt_dosurf )
+
+          ntdim_surf(1) = dosurf_time_count(1) + CEILING(                      &
+                        ( end_time - MAX(                                      &
+                            MERGE(   skip_time_dosurf_av, skip_time_dosurf_av  &
+                                  + spinup_time, data_output_during_spinup ),  &
+                            simulated_time_at_begin )                          &
+                        ) / dt_dosurf_av )
+
+!
+!--       Create NetCDF4 files for parallel writing
+          DO  av = 0, 1
+!
+!--          If there is no instantaneous data (av=0) or averaged data (av=1)
+!--          requested, do not create the corresponding NetCDF file
+             IF ( dosurf_no(av) == 0 ) CYCLE
+
+             IF ( av == 0 )  THEN
+                filename = 'SURFACE_DATA_NETCDF' // TRIM( coupling_char )
+             ELSE
+                filename = 'SURFACE_DATA_AV_NETCDF' // TRIM( coupling_char )
+             ENDIF
+!
+!--          Open file using netCDF4/HDF5 format, parallel
+             nc_stat = NF90_CREATE( TRIM(filename),                            &
+                                    IOR( NF90_NOCLOBBER,                       &
+                                    IOR( NF90_NETCDF4, NF90_MPIIO ) ),         &
+                                    id_set_surf(av), COMM = comm2d, INFO = MPI_INFO_NULL )
+             CALL netcdf_handle_error( 'surface_data_output_mod', 5550 )
+
+             !- Write some global attributes
+             IF ( av == 0 )  THEN
+                CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data', TRIM( run_description_header ), 5551 )
+                time_average_text = ' '
+             ELSE
+                CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data_av', TRIM( run_description_header ), 5552 )
+                WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_surf
+                nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'time_avg',  &
+                                        TRIM( time_average_text ) )
+                CALL netcdf_handle_error( 'surface_data_output_mod', 5553 )
+             ENDIF
+
+
+!
+!--          Define time coordinate for surface data.
+!--          For parallel output the time dimension has to be limited (ntdim_surf),
+!--          otherwise the performance drops significantly.
+             CALL netcdf_create_dim( id_set_surf(av), 'time', ntdim_surf(av),  &
+                                     id_dim_time_surf(av), 5554 )
+
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_time_surf(av) /),    &
+                                     'time', NF90_DOUBLE, id_var_time_surf(av),      &
+                                     'seconds since '//TRIM(init_model%origin_time), &
+                                     'time', 5555, 5555, 5555 )
+             CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'standard_name', 'time', 5556)
+             CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'axis', 'T', 5557)
+!
+!--          Define spatial dimensions and coordinates:
+!--          Define index of surface element
+             CALL netcdf_create_dim( id_set_surf(av), 's', surfaces%ns_total,  &
+                                     id_dim_s_surf(av), 5558 )
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
+                                     's', NF90_DOUBLE, id_var_s_surf(av),      &
+                                     '1', '', 5559, 5559, 5559 )
+!
+!--          Define x coordinate
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
+                                     'xs', NF90_DOUBLE, id_var_xs_surf(av),    &
+                                     'meters', '', 5561, 5561, 5561 )
+!
+!--           Define y coordinate
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
+                                     'ys', NF90_DOUBLE, id_var_ys_surf(av),    &
+                                     'meters', '', 5562, 5562, 5562 )
+!
+!--          Define z coordinate
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
+                                     'zs', NF90_DOUBLE, id_var_zs_surf(av),    &
+                                     'meters', '', 5560, 5560, 5560 )
+
+!
+!--          Define UTM coordinates
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /),    &
+                                     'Es_UTM', NF90_DOUBLE, id_var_etum_surf(av), &
+                                     'meters', '', 5563, 5563, 5563 )
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /),    &
+                                     'Ns_UTM', NF90_DOUBLE, id_var_nutm_surf(av), &
+                                     'meters', '', 5564, 5564, 5564 )
+!
+!--          Define the variables
+             var_list = ';'
+             i = 1
+
+             DO WHILE ( dosurf(av,i)(1:1) /= ' ' )
+
+                CALL netcdf_create_var( id_set_surf(av),(/ id_dim_s_surf(av),  &
+                                        id_dim_time_surf(av) /), dosurf(av,i), &
+                                        NF90_REAL4, id_var_dosurf(av,i),       &
+                                        dosurf_unit(av,i), dosurf(av,i), 5565, &
+                                        5565, 5565, .TRUE. )
+!
+!--                Set no fill for every variable to increase performance.
+                nc_stat = NF90_DEF_VAR_FILL( id_set_surf(av),     &
+                                             id_var_dosurf(av,i), &
+                                             1, 0 )
+                CALL netcdf_handle_error( 'surface_data_output_init', 5566 )
+!
+!--                Set collective io operations for parallel io
+                nc_stat = NF90_VAR_PAR_ACCESS( id_set_surf(av),     &
+                                               id_var_dosurf(av,i), &
+                                               NF90_COLLECTIVE )
+                CALL netcdf_handle_error( 'surface_data_output_init', 5567 )
+                var_list = TRIM( var_list ) // TRIM( dosurf(av,i) ) // ';'
+
+                i = i + 1
+
+             ENDDO
+!
+!--          Write the list of variables as global attribute (this is used by
+!--          restart runs and by combine_plot_fields)
+             nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'VAR_LIST', &
+                                     var_list )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5568 )
+
+!
+!--          Set general no fill, otherwise the performance drops significantly for
+!--          parallel output.
+             nc_stat = NF90_SET_FILL( id_set_surf(av), NF90_NOFILL, oldmode )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5569 )
+
+!
+!--          Leave netCDF define mode
+             nc_stat = NF90_ENDDEF( id_set_surf(av) )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5570 )
+
+!
+!--          These data are only written by PE0 for parallel output to increase
+!--          the performance.
+             IF ( myid == 0 )  THEN
+!
+!--             Write data for surface indices
+                ALLOCATE( netcdf_data_1d(1:surfaces%ns_total) )
+
+                DO  i = 1, surfaces%ns_total
+                   netcdf_data_1d(i) = i
+                ENDDO
+
+                nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_s_surf(av),  &
+                                        netcdf_data_1d, start = (/ 1 /),     &
+                                        count = (/ surfaces%ns_total /) )
+                CALL netcdf_handle_error( 'surface_data_output_init', 5571 )
+
+                DEALLOCATE( netcdf_data_1d )
+
+             ENDIF
+
+!
+!--          Write surface positions to file
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_xs_surf(av),      &
+                                     surfaces%xs, start = (/ surfaces%s(1) /), &
+                                     count = (/ surfaces%ns /) )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5572 )
+
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_ys_surf(av),      &
+                                     surfaces%ys, start = (/ surfaces%s(1) /), &
+                                     count = (/ surfaces%ns /) )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5573 )
+
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_zs_surf(av),      &
+                                     surfaces%zs, start = (/ surfaces%s(1) /), &
+                                     count = (/ surfaces%ns /) )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5574 )
+
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_etum_surf(av),        &
+                                     surfaces%es_utm, start = (/ surfaces%s(1) /), &
+                                     count = (/ surfaces%ns /) )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5575 )
+
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_nutm_surf(av),        &
+                                     surfaces%ns_utm, start = (/ surfaces%s(1) /), &
+                                     count = (/ surfaces%ns /) )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5576 )
+             
+          ENDDO
+#endif
+
+       ENDIF
 
    END SUBROUTINE surface_data_output_init
@@ -892,25 +1299,55 @@
       IF ( dosurf_no(av) == 0 )  RETURN
 !
-!--   Open files
-      CALL check_open( 25+av )
-!
-!--   Write coordinates 
-      IF ( .NOT. first_output(av) )  THEN
-         DO  i = 0, io_blocks-1
-            IF ( i == io_group )  THEN
-               WRITE ( 25+av )  surfaces%npoints
-               WRITE ( 25+av )  surfaces%npoints_total
-               WRITE ( 25+av )  surfaces%ns
-               WRITE ( 25+av )  surfaces%ns_total
-               WRITE ( 25+av )  surfaces%points
-               WRITE ( 25+av )  surfaces%polygons
-            ENDIF
+!--   In case of VTK output, check if binary files are open and write coordinates. 
+      IF ( to_vtk )  THEN
+
+         CALL check_open( 25+av )
+
+         IF ( .NOT. first_output(av) )  THEN
+            DO  i = 0, io_blocks-1
+               IF ( i == io_group )  THEN
+                  WRITE ( 25+av )  surfaces%npoints
+                  WRITE ( 25+av )  surfaces%npoints_total
+                  WRITE ( 25+av )  surfaces%ns
+                  WRITE ( 25+av )  surfaces%ns_total
+                  WRITE ( 25+av )  surfaces%points
+                  WRITE ( 25+av )  surfaces%polygons
+               ENDIF
 #if defined( __parallel )
-            CALL MPI_BARRIER( comm2d, ierr )
+               CALL MPI_BARRIER( comm2d, ierr )
 #endif
-            first_output(av) = .TRUE.
-         ENDDO
+               first_output(av) = .TRUE.
+            ENDDO
+         ENDIF
       ENDIF
-
+!
+!--   In case of NetCDF output, check if enough time steps are available in file
+!--   and update time axis.
+      IF ( to_netcdf )  THEN
+         IF ( dosurf_time_count(av) + 1 > ntdim_surf(av) )  THEN
+            WRITE ( message_string, * ) 'Output of surface data is not given at t=',          &
+                                        time_since_reference_point, 's because the maximum ', & 
+                                        'number of output time levels is exceeded.'
+            CALL message( 'surface_data_output', 'PA0???', 0, 1, 0, 6, 0 )
+            
+            RETURN
+
+         ENDIF
+!
+!--      Update the netCDF time axis
+!--      In case of parallel output, this is only done by PE0 to increase the
+!--      performance.
+         dosurf_time_count(av) = dosurf_time_count(av) + 1
+         IF ( myid == 0 )  THEN
+            nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_time_surf(av),  &
+                                    (/ time_since_reference_point /),       &
+                                    start = (/ dosurf_time_count(av) /),    &
+                                    count = (/ 1 /) )
+            CALL netcdf_handle_error( 'surface_data_output', 6666 )
+         ENDIF
+      ENDIF
+
+!
+!--   Cycle through output quantities and write them to file.
       n_out = 0
       DO  WHILE ( dosurf(av,n_out+1)(1:1) /= ' ' )
@@ -2269,17 +2706,30 @@
 !--      - Dimension: 1-nsurfaces, 1-npoints - can be ordered consecutively
 !--      - Distinguish between average and non-average data
-         DO  i = 0, io_blocks-1
-            IF ( i == io_group )  THEN
-               WRITE ( 25+av )  LEN_TRIM( 'time' )
-               WRITE ( 25+av )  'time'
-               WRITE ( 25+av )  time_since_reference_point
-               WRITE ( 25+av )  LEN_TRIM( trimvar )
-               WRITE ( 25+av )  TRIM( trimvar )
-               WRITE ( 25+av )  surfaces%var_out
-            ENDIF
+         IF ( to_vtk )  THEN
+            DO  i = 0, io_blocks-1
+               IF ( i == io_group )  THEN
+                  WRITE ( 25+av )  LEN_TRIM( 'time' )
+                  WRITE ( 25+av )  'time'
+                  WRITE ( 25+av )  time_since_reference_point
+                  WRITE ( 25+av )  LEN_TRIM( trimvar )
+                  WRITE ( 25+av )  TRIM( trimvar )
+                  WRITE ( 25+av )  surfaces%var_out
+               ENDIF
 #if defined( __parallel )
-            CALL MPI_BARRIER( comm2d, ierr )
+               CALL MPI_BARRIER( comm2d, ierr )
 #endif
-         ENDDO
+            ENDDO
+         ENDIF
+         
+         IF ( to_netcdf )  THEN
+!
+!--         Write output array to file
+            nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_dosurf(av,n_out),            &
+                                    surfaces%var_out,                                    &
+                                    start = (/ surfaces%s(1), dosurf_time_count(av) /),  &
+                                    count = (/ surfaces%ns, 1 /) )
+            CALL netcdf_handle_error( 'surface_data_output', 6667 )
+         
+         ENDIF
 
       ENDDO
@@ -3553,5 +4003,6 @@
                                   averaging_interval_surf, data_output_surf,   &
                                   dt_dosurf, dt_dosurf_av,                     &
-                                  skip_time_dosurf, skip_time_dosurf_av
+                                  skip_time_dosurf, skip_time_dosurf_av,       &
+                                  to_netcdf, to_vtk
                                   
        line = ' '
@@ -3599,4 +4050,5 @@
 
        CHARACTER(LEN=100) ::  trimvar !< dummy for single output variable
+       CHARACTER(LEN=100) ::  unit    !< dummy for unit of output variable
 
        INTEGER(iwp) ::  av    !< id indicating average or non-average data output
@@ -3621,5 +4073,14 @@
                         'PA0087', 1, 2, 0, 6, 0 )
        ENDIF
-
+       
+#if ! defined( __netcdf4_parallel )
+!
+!--    Surface output via NetCDF requires parallel NetCDF
+       IF ( to_netcdf )  THEN
+          message_string = 'to_netcdf = .True. requires parallel NetCDF'
+          CALL message( 'surface_data_output_check_parameters',                &
+                        'PA0116', 1, 2, 0, 6, 0 )
+       ENDIF
+#endif
 !
 !--   Count number of output variables and separate output strings for 
@@ -3646,5 +4107,6 @@
 
 !
-!--      Check if all output variables are known
+!--      Check if all output variables are known and assign a unit
+         unit = 'not set'
          SELECT CASE ( TRIM( trimvar ) )
 
@@ -3655,31 +4117,49 @@
                              'PA0087', 1, 2, 0, 6, 0 )
 
-            CASE ( 'us', 'ts', 'qs', 'ss', 'qcs', 'ncs', 'qrs', 'nrs', 'ol' )
-
-            CASE ( 'z0', 'z0h', 'z0q' )
-
-            CASE ( 'theta1', 'qv1', 'thetav1' )
+            CASE ( 'us', 'uvw1' )
+               unit = 'm/s'
+
+            CASE ( 'ss', 'qcs', 'ncs', 'qrs', 'nrs' )
+               unit = '1'
+
+            CASE ( 'z0', 'z0h', 'z0q', 'ol' )
+               unit = 'm'
+
+            CASE ( 'ts', 'theta1', 'thetav1', 'theta_surface', 'thetav_surface' )
+               unit = 'K'
 
             CASE ( 'usws', 'vsws' )
+               unit = 'm2/s2'
             
-            CASE ( 'uvw1' )
-
             CASE ( 'qcsws', 'ncsws', 'qrsws', 'nrsws', 'sasws' )
 
-            CASE ( 'shf', 'qsws', 'ssws' )
-
-            CASE ( 'q_surface', 'theta_surface', 'thetav_surface' )
+            CASE ( 'shf' )
+               unit = 'K m/s'
+
+            CASE ( 'qsws' )
+               unit = 'kg/kg m/s'
+
+            CASE ( 'ssws' )
+               unit = 'kg/m2/s'
+
+            CASE ( 'qs', 'q_surface', 'qv1' )
+               unit = 'kg/kg'
 
             CASE ( 'rad_net' )
+               unit = 'W/m2'
 
             CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_dif', 'rad_lw_ref',     &
                    'rad_lw_res' )
+               unit = 'W/m2'
 
             CASE ( 'rad_sw_in', 'rad_sw_out', 'rad_sw_dif', 'rad_sw_ref',     &
                    'rad_sw_res', 'rad_sw_dir' )
+               unit = 'W/m2'
 
             CASE ( 'ghf' )
+               unit = 'W/m2'
 
             CASE ( 'r_a', 'r_canopy', 'r_soil', 'r_s' )
+               unit = 's/m'
 
             CASE DEFAULT
@@ -3689,4 +4169,6 @@
                              'PA0087', 1, 2, 0, 6, 0 )
          END SELECT
+
+         dosurf_unit(av,dosurf_no(av)) = unit
 
        ENDDO
@@ -3717,17 +4199,20 @@
       IF ( dosurf_no(av) == 0 )  RETURN
 !
-!--   Open files
-      CALL check_open( 25+av )
-!
-!--   Write time coordinate
-      DO  i = 0, io_blocks-1
-         IF ( i == io_group )  THEN
-            WRITE ( 25+av )  LEN_TRIM( 'END' )
-            WRITE ( 25+av )  'END'
-         ENDIF
+!--   If output to VTK files is enabled, check if files are open and write
+!--   an end-of-file statement.
+      IF ( to_vtk )  THEN
+         CALL check_open( 25+av )
+!
+!--      Write time coordinate
+         DO  i = 0, io_blocks-1
+            IF ( i == io_group )  THEN
+               WRITE ( 25+av )  LEN_TRIM( 'END' )
+               WRITE ( 25+av )  'END'
+            ENDIF
 #if defined( __parallel )
-         CALL MPI_BARRIER( comm2d, ierr )
+            CALL MPI_BARRIER( comm2d, ierr )
 #endif
-      ENDDO    
+         ENDDO    
+      ENDIF
 
     END SUBROUTINE surface_data_output_last_action    
