Index: /palm/trunk/SCRIPTS/.palm.iofiles
===================================================================
--- /palm/trunk/SCRIPTS/.palm.iofiles	(revision 3420)
+++ /palm/trunk/SCRIPTS/.palm.iofiles	(revision 3421)
@@ -52,4 +52,7 @@
 LIST_PROFIL*               out:tr   *         $base_data/$run_identifier/MONITORING     _list_pr
 #
+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
+#
 DATA_1D_PR_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _pr        nc
 DATA_1D_SP_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _sp        nc
Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 3420)
+++ /palm/trunk/SOURCE/Makefile	(revision 3421)
@@ -25,4 +25,9 @@
 # -----------------
 # $Id$
+# 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
 # 
@@ -611,4 +616,5 @@
 	surface_layer_fluxes_mod.f90 \
 	surface_mod.f90 \
+	surface_output_mod.f90 \
 	swap_timelevel.f90 \
 	synthetic_turbulence_generator_mod.f90 \
@@ -798,4 +804,5 @@
 	spectra_mod.o \
 	subsidence_mod.o \
+	surface_output_mod.o \
 	synthetic_turbulence_generator_mod.o \
 	turbulence_closure_mod.o \
@@ -1052,4 +1059,5 @@
 	surface_layer_fluxes_mod.o \
 	surface_mod.o \
+	surface_output_mod.o \
 	synthetic_turbulence_generator_mod.o \
 	turbulence_closure_mod.o \
@@ -1289,4 +1297,5 @@
 	mod_kinds.o \
 	modules.o \
+	netcdf_data_input_mod.o \
 	ocean_mod.o \
 	radiation_model_mod.o \
@@ -1358,4 +1367,5 @@
 	read_restart_data_mod.o \
 	spectra_mod.o \
+	surface_output_mod.o \
 	synthetic_turbulence_generator_mod.o \
 	turbulence_closure_mod.o \
@@ -1574,4 +1584,8 @@
 	mod_kinds.o \
 	modules.o
+surface_output_mod.o: \
+	mod_kinds.o \
+	modules.o \
+	surface_mod.o
 swap_timelevel.o: \
 	bulk_cloud_model_mod.o \
@@ -1625,4 +1639,5 @@
 	surface_layer_fluxes_mod.o \
 	surface_mod.o \
+	surface_output_mod.o \
 	synthetic_turbulence_generator_mod.o \
 	time_to_string.o \
Index: /palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/average_3d_data.f90	(revision 3420)
+++ /palm/trunk/SOURCE/average_3d_data.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3419 2018-10-24 17:27:31Z gronemeier
 ! (from resler branch)
 ! Add biometeorology output,
@@ -250,5 +253,5 @@
              ENDIF
 
-          CASE ( 'lpt' )
+          CASE ( 'thetal' )
              IF ( ALLOCATED( lpt_av ) ) THEN
                 DO  i = nxlg, nxrg
@@ -313,5 +316,5 @@
              ENDIF
 
-          CASE ( 'pt' )
+          CASE ( 'theta' )
              IF ( ALLOCATED( pt_av ) ) THEN
                 DO  i = nxlg, nxrg
@@ -463,5 +466,5 @@
              ENDIF
 
-          CASE ( 'u*' )
+          CASE ( 'us*' )
              IF ( ALLOCATED( us_av ) ) THEN
                 DO  i = nxlg, nxrg
@@ -484,5 +487,5 @@
              ENDIF
 
-          CASE ( 'vpt' )
+          CASE ( 'thetav' )
              IF ( ALLOCATED( vpt_av ) ) THEN
                 DO  i = nxlg, nxrg
Index: /palm/trunk/SOURCE/check_open.f90
===================================================================
--- /palm/trunk/SOURCE/check_open.f90	(revision 3420)
+++ /palm/trunk/SOURCE/check_open.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Open files for surface data
+! 
+! 3399 2018-10-22 20:13:45Z gronemeier
 ! unused variables and format statements removed
 ! 
@@ -455,4 +458,17 @@
 
           ENDIF
+          
+       CASE ( 25 )
+!
+!--       Binary files for surface data
+          OPEN ( 25, FILE='SURFACE_DATA_BIN'//TRIM( coupling_char )//myid_char,  &
+                 FORM='UNFORMATTED', POSITION='APPEND' )
+
+       CASE ( 26 )
+!
+!--       Binary files for averaged surface data
+          OPEN ( 26, FILE='SURFACE_DATA_AV_BIN'//TRIM( coupling_char )//myid_char,  &
+                 FORM='UNFORMATTED', POSITION='APPEND' )
+       
 
        CASE ( 30 )
@@ -1131,4 +1147,5 @@
           ENDIF
 
+
 !
 !--    nc-file for virtual flight measurements
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 3420)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 3421)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! Add checks for surface data output
+! 
+! 3419 2018-10-24 17:27:31Z gronemeier
 ! Remove offline nesting in if clause for pressure top boundary condition
 ! 
@@ -777,4 +781,7 @@
 
     USE statistics
+
+    USE surface_output_mod,                                                    &
+        ONLY:  surface_output_check_parameters
 
     USE synthetic_turbulence_generator_mod,                                    &
@@ -1494,4 +1501,7 @@
     IF ( urban_surface )        CALL usm_check_parameters
     IF ( wind_turbine )         CALL wtm_check_parameters
+!
+!-- Checks for surface data output
+    IF ( surface_data_output )  CALL surface_output_check_parameters
 
 !
@@ -2363,5 +2373,5 @@
              hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'pt', '#pt' )
+          CASE ( 'theta', '#theta' )
              IF ( .NOT. bulk_cloud_model ) THEN
                 dopr_index(i) = 4
@@ -2452,15 +2462,15 @@
              hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'w"pt"' )
+          CASE ( 'w"theta"' )
              dopr_index(i) = 16
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'w*pt*' )
+          CASE ( 'w*theta*' )
              dopr_index(i) = 17
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'wpt' )
+          CASE ( 'wtheta' )
              dopr_index(i) = 18
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
@@ -2479,10 +2489,10 @@
              IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
 
-          CASE ( 'w*pt*BC' )
+          CASE ( 'w*theta*BC' )
              dopr_index(i) = 21
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'wptBC' )
+          CASE ( 'wthetaBC' )
              dopr_index(i) = 22
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
@@ -2504,5 +2514,5 @@
              hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'pt*2' )
+          CASE ( 'theta*2' )
              dopr_index(i) = 33
              dopr_unit(i)  = 'K2'
@@ -2603,5 +2613,5 @@
              ENDIF
 
-          CASE ( 'lpt', '#lpt' )
+          CASE ( 'thetal', '#thetal' )
              IF ( .NOT. bulk_cloud_model ) THEN
                 message_string = 'data_output_pr = ' //                        &
@@ -2621,5 +2631,5 @@
              ENDIF
 
-          CASE ( 'vpt', '#vpt' )
+          CASE ( 'thetav', '#thetav' )
              dopr_index(i) = 44
              dopr_unit(i)  = 'K'
@@ -2632,15 +2642,15 @@
              ENDIF
 
-          CASE ( 'w"vpt"' )
+          CASE ( 'w"thetav"' )
              dopr_index(i) = 45
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'w*vpt*' )
+          CASE ( 'w*thetav*' )
              dopr_index(i) = 46
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'wvpt' )
+          CASE ( 'wthetav' )
              dopr_index(i) = 47
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
@@ -2798,30 +2808,30 @@
              hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'u"pt"' )
+          CASE ( 'u"theta"' )
              dopr_index(i) = 58
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'u*pt*' )
+          CASE ( 'u*theta*' )
              dopr_index(i) = 59
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'upt_t' )
+          CASE ( 'utheta_t' )
              dopr_index(i) = 60
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'v"pt"' )
+          CASE ( 'v"theta"' )
              dopr_index(i) = 61
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'v*pt*' )
+          CASE ( 'v*theta*' )
              dopr_index(i) = 62
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'vpt_t' )
+          CASE ( 'thetav_t' )
              dopr_index(i) = 63
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
@@ -2855,10 +2865,10 @@
              hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'rho_air' )
+          CASE ( 'rho' )
              dopr_index(i)  = 119
              dopr_unit(i)   = 'kg/m3'
              hom(:,2,119,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'rho_air_zw' )
+          CASE ( 'rho_zw' )
              dopr_index(i)  = 120
              dopr_unit(i)   = 'kg/m3'
@@ -3023,5 +3033,5 @@
              unit = 'm2/s2'
 
-          CASE ( 'lpt' )
+          CASE ( 'thetal' )
              IF (  .NOT.  bulk_cloud_model )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3041,5 +3051,5 @@
              IF ( TRIM( var ) == 'pr' )  unit = 'm'
 
-          CASE ( 'q', 'vpt' )
+          CASE ( 'q', 'thetav' )
              IF (  .NOT.  humidity )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3048,5 +3058,5 @@
              ENDIF
              IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
-             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
+             IF ( TRIM( var ) == 'thetav' )  unit = 'K'
 
           CASE ( 'ql' )
@@ -3084,7 +3094,7 @@
              unit = 'kg/m3'
 
-          CASE ( 'p', 'pt', 'u', 'v', 'w' )
+          CASE ( 'p', 'theta', 'u', 'v', 'w' )
              IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
-             IF ( TRIM( var ) == 'pt' )  unit = 'K'
+             IF ( TRIM( var ) == 'theta' )  unit = 'K'
              IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
              IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
@@ -3093,5 +3103,5 @@
 
           CASE ( 'ghf*', 'lwp*', 'ol*', 'qsws*', 'r_a*',                       &
-                 'shf*', 'ssws*', 't*', 'tsurf*', 'u*', 'z0*', 'z0h*', 'z0q*' )
+                 'shf*', 'ssws*', 't*', 'tsurf*', 'us*', 'z0*', 'z0h*', 'z0q*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                 message_string = 'illegal value for data_output: "' //         &
@@ -3141,5 +3151,5 @@
              IF ( TRIM( var ) == 't*'     )  unit = 'K'
              IF ( TRIM( var ) == 'tsurf*' )  unit = 'K'  
-             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
+             IF ( TRIM( var ) == 'us*'    )  unit = 'm/s'
              IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
              IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
Index: /palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_2d.f90	(revision 3420)
+++ /palm/trunk/SOURCE/data_output_2d.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3419 2018-10-24 17:27:31Z gronemeier
 ! minor formatting (kanani)
 ! chem_data_output_2d subroutine added (basit)
@@ -562,5 +565,5 @@
                 IF ( mode == 'xy' )  level_z = zu
 
-             CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
+             CASE ( 'thetal_xy', 'thetal_xz', 'thetal_yz' )
                 IF ( av == 0 )  THEN
                    to_be_resorted => pt
@@ -750,5 +753,5 @@
                 ENDIF
 
-             CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
+             CASE ( 'theta_xy', 'theta_xz', 'theta_yz' )
                 IF ( av == 0 )  THEN
                    IF ( .NOT. bulk_cloud_model ) THEN
@@ -1121,5 +1124,5 @@
                 ENDIF
 
-             CASE ( 'u*_xy' )        ! 2d-array
+             CASE ( 'us*_xy' )        ! 2d-array
                 IF ( av == 0 )  THEN
                    DO  m = 1, surf_def_h(0)%ns
@@ -1172,5 +1175,5 @@
                 ENDIF
 
-             CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
+             CASE ( 'thetav_xy', 'thetav_xz', 'thetav_yz' )
                 IF ( av == 0 )  THEN
                    to_be_resorted => vpt
Index: /palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_3d.f90	(revision 3420)
+++ /palm/trunk/SOURCE/data_output_3d.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3419 2018-10-24 17:27:31Z gronemeier
 ! bugfix: nx, ny are required in non-parallel case
 ! 
@@ -418,5 +421,5 @@
              ENDIF
 
-          CASE ( 'lpt' )
+          CASE ( 'thetal' )
              IF ( av == 0 )  THEN
                 to_be_resorted => pt
@@ -512,5 +515,5 @@
              ENDIF
 
-          CASE ( 'pt' )
+          CASE ( 'theta' )
              IF ( av == 0 )  THEN
                 IF ( .NOT. bulk_cloud_model ) THEN
@@ -671,5 +674,5 @@
              ENDIF
 
-          CASE ( 'vpt' )
+          CASE ( 'thetav' )
              IF ( av == 0 )  THEN
                 to_be_resorted => vpt
Index: /palm/trunk/SOURCE/data_output_dvrp.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_dvrp.f90	(revision 3420)
+++ /palm/trunk/SOURCE/data_output_dvrp.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3419 2018-10-24 17:27:31Z gronemeier
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -337,5 +340,5 @@
                    ENDDO
 
-                CASE ( 'pt', 'pt_xy', 'pt_xz', 'pt_yz' )
+                CASE ( 'theta', 'theta_xy', 'theta_xz', 'theta_yz' )
                    IF ( .NOT. bulk_cloud_model ) THEN
                       DO  i = nxl_dvrp, nxr_dvrp+1
@@ -409,5 +412,5 @@
                    ENDIF
 
-                CASE ( 'u*_xy' )
+                CASE ( 'us*_xy' )
                    DO  i = nxl_dvrp, nxr_dvrp+1
                       DO  j = nys_dvrp, nyn_dvrp+1
Index: /palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_mask.f90	(revision 3420)
+++ /palm/trunk/SOURCE/data_output_mask.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3419 2018-10-24 17:27:31Z gronemeier
 ! Minor formatting (kanani)
 ! Call for chem_data_output_mask added (basit)
@@ -273,5 +276,5 @@
              ENDIF
 
-          CASE ( 'lpt' )
+          CASE ( 'thetal' )
              IF ( av == 0 )  THEN
                 to_be_resorted => pt
@@ -365,5 +368,5 @@
              ENDIF
 
-          CASE ( 'pt' )
+          CASE ( 'theta' )
              IF ( av == 0 )  THEN
                 IF ( .NOT. bulk_cloud_model ) THEN
@@ -481,5 +484,5 @@
              ENDIF
 
-          CASE ( 'rho_ocean' )
+          CASE ( 'rho_sea_water' )
              IF ( av == 0 )  THEN
                 to_be_resorted => rho_ocean
@@ -516,5 +519,5 @@
              ENDIF
 
-          CASE ( 'vpt' )
+          CASE ( 'thetav' )
              IF ( av == 0 )  THEN
                 to_be_resorted => vpt
@@ -531,4 +534,5 @@
 
           CASE DEFAULT
+
 !
 !--          Radiation quantity
Index: /palm/trunk/SOURCE/data_output_spectra.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_spectra.f90	(revision 3420)
+++ /palm/trunk/SOURCE/data_output_spectra.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3419 2018-10-24 17:27:31Z gronemeier
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -179,5 +182,5 @@
                 pr = 3
 
-             CASE ( 'pt' )
+             CASE ( 'theta' )
                 pr = 4
 
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 3420)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Initialize surface data output
+! 
+! 3415 2018-10-24 11:57:50Z suehring
 ! Set bottom boundary condition for geostrophic wind components in inifor
 ! initialization
@@ -625,4 +628,7 @@
         ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
                 surf_usm_h, get_topography_top_index_ji, vertical_surfaces_exist
+                
+    USE surface_output_mod,                                                    &
+        ONLY:  surface_output_init
    
     USE transpose_indices
@@ -2400,4 +2406,9 @@
        CALL gust_init( dots_label, dots_unit, dots_num, dots_max )
     ENDIF
+!
+!-- Initialize surface data output
+    IF ( surface_data_output )  THEN
+       CALL surface_output_init
+    ENDIF
 
 !
Index: /palm/trunk/SOURCE/init_masks.f90
===================================================================
--- /palm/trunk/SOURCE/init_masks.f90	(revision 3420)
+++ /palm/trunk/SOURCE/init_masks.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3419 2018-10-24 17:27:31Z gronemeier
 ! ocean renamed ocean_mode
 ! 
@@ -316,5 +319,5 @@
                 unit = 'm2/s2'
 
-             CASE ( 'lpt' )
+             CASE ( 'thetal' )
                 IF ( .NOT. bulk_cloud_model )  THEN
                    WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
@@ -358,5 +361,5 @@
                 IF ( TRIM( var ) == 'pr' )  unit = 'm'
 
-             CASE ( 'q', 'vpt' )
+             CASE ( 'q', 'thetav' )
                 IF ( .NOT. humidity )  THEN
                    WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
@@ -365,5 +368,5 @@
                 ENDIF
                 IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
-                IF ( TRIM( var ) == 'vpt' )  unit = 'K'
+                IF ( TRIM( var ) == 'thetav' )  unit = 'K'
 
              CASE ( 'qc' )
@@ -414,5 +417,5 @@
                 unit = 'kg/kg'
 
-             CASE ( 'rho_ocean' )
+             CASE ( 'rho_sea_water' )
                 IF ( .NOT. ocean_mode )  THEN
                    WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
@@ -438,5 +441,5 @@
                 unit = 'psu'
 
-             CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*', 'z0h*' )
+             CASE ( 'us*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*', 'z0h*' )
                 WRITE ( message_string, * ) 'illegal value for data_',         &
                      'output: "', TRIM( var ), '" is only allowed',            &
@@ -444,7 +447,7 @@
                 CALL message( 'init_masks', 'PA0111', 1, 2, 0, 6, 0 )
 
-             CASE ( 'p', 'pt', 'u', 'v', 'w' )
+             CASE ( 'p', 'theta', 'u', 'v', 'w' )
                 IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
-                IF ( TRIM( var ) == 'pt' )  unit = 'K'
+                IF ( TRIM( var ) == 'theta' )  unit = 'K'
                 IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
                 IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 3420)
+++ /palm/trunk/SOURCE/modules.f90	(revision 3421)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! +surface_data_output
+! 
+! 3355 2018-10-16 14:03:34Z knoop
 ! (from branch resler)
 ! Increase dimension of uv_heights etc.
@@ -1353,4 +1357,5 @@
     LOGICAL ::  sloping_surface = .FALSE.                        !< use sloped surface? (namelist parameter alpha_surface)
     LOGICAL ::  spinup = .FALSE.                                 !< perform model spinup without atmosphere code?
+    LOGICAL ::  surface_data_output = .FALSE.                    !< output of surface data
     LOGICAL ::  stop_dt = .FALSE.                                !< internal switch to stop the time stepping
     LOGICAL ::  synchronous_exchange = .FALSE.                   !< namelist parameter
@@ -2027,20 +2032,25 @@
     INTEGER(iwp), PARAMETER ::  crmax = 100  !< maximum number of coordinate systems for profile output
 
-    CHARACTER (LEN=20), DIMENSION(20) ::  cross_ts_profiles = &  !< time series to be plotted into one coordinate system, respectively
-                           (/ ' E E*               ', ' dt                 ', &
-                              ' u* w*              ', ' th*                ', &
-                              ' umax vmax wmax     ', ' div_old div_new    ', &
-                              ' z_i_wpt z_i_pt     ', ' w"pt"0 w"pt" wpt   ', &
-                              ' pt(0) pt(zp)       ', ' splux spluy spluz  ', &
-                              ' L                  ',                         &
-                            ( '                    ', i9 = 1, 9 ) /)
+    CHARACTER (LEN=27), DIMENSION(20) ::  cross_ts_profiles = &  !< time series to be plotted into one coordinate system, respectively
+                           (/ ' E E*                      ', &
+                              ' dt                        ', &
+                              ' u* w*                     ', &
+                              ' th*                       ', &
+                              ' umax vmax wmax            ', &
+                              ' div_old div_new           ', &
+                              ' zi_wtheta zi_theta        ', &
+                              ' w"theta"0 w"theta" wtheta ', &
+                              ' theta(0) theta(zp)        ', &
+                              ' splux spluy spluz         ', &
+                              ' L                         ', &
+                            ( '                           ', i9 = 1, 9 ) /)
 
     CHARACTER (LEN=100), DIMENSION(crmax) ::  cross_profiles = &  !< quantities to be plotted into one coordinate system, respectively
-                           (/ ' u v                           ', &
-                              ' pt                            ', &
-                              ' w"pt" w*pt* w*pt*BC wpt wptBC ', &
-                              ' w"u" w*u* wu w"v" w*v* wv     ', &
-                              ' km kh                         ', &
-                              ' l                             ', &
+                          (/ ' u v                                          ', &
+                             ' pt                                           ', &
+                             ' w"theta" w*theta* w*theta*BC wtheta wthetaBC ', &
+                             ' w"u" w*u* wu w"v" w*v* wv                    ', &
+                             ' km kh                                        ', &
+                             ' l                                            ', &
                          ( '                               ', i9 = 1, 94 ) /)
 
Index: /palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3420)
+++ /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! add default values of origin_x/y/z
+! 
+! 3404 2018-10-23 13:29:11Z suehring
 ! Consider time-dependent geostrophic wind components in offline nesting
 ! 
@@ -328,17 +331,17 @@
        LOGICAL ::  from_file_w      = .FALSE. !< flag indicating whether w is already initialized from file
 
-       REAL(wp) ::  fill_msoil       !< fill value for soil moisture
-       REAL(wp) ::  fill_pt          !< fill value for pt
-       REAL(wp) ::  fill_q           !< fill value for q
-       REAL(wp) ::  fill_tsoil       !< fill value for soil temperature
-       REAL(wp) ::  fill_u           !< fill value for u
-       REAL(wp) ::  fill_v           !< fill value for v
-       REAL(wp) ::  fill_w           !< fill value for w
-       REAL(wp) ::  latitude         !< latitude of the southern model boundary
-       REAL(wp) ::  longitude        !< longitude of the western model boundary
-       REAL(wp) ::  origin_x         !< x position of the western model boundary
-       REAL(wp) ::  origin_y         !< y position of the northern model boundary
-       REAL(wp) ::  origin_z         !< reference height of input data
-       REAL(wp) ::  rotation_angle   !< rotation angle of input data
+       REAL(wp) ::  fill_msoil              !< fill value for soil moisture
+       REAL(wp) ::  fill_pt                 !< fill value for pt
+       REAL(wp) ::  fill_q                  !< fill value for q
+       REAL(wp) ::  fill_tsoil              !< fill value for soil temperature
+       REAL(wp) ::  fill_u                  !< fill value for u
+       REAL(wp) ::  fill_v                  !< fill value for v
+       REAL(wp) ::  fill_w                  !< fill value for w
+       REAL(wp) ::  latitude                !< latitude of the southern model boundary
+       REAL(wp) ::  longitude               !< longitude of the western model boundary
+       REAL(wp) ::  origin_x = 0.0_wp       !< x position of the western model boundary (default: Hannover)
+       REAL(wp) ::  origin_y = 0.0_wp       !< y position of the northern model boundary (default: Hannover)
+       REAL(wp) ::  origin_z = 0.0_wp       !< reference height of input data (default: Hannover)
+       REAL(wp) ::  rotation_angle = 0.0_wp !< rotation angle of input data
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  msoil_1d     !< initial vertical profile of soil moisture
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3420)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3421)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Bugfix: move ocean output variables to ocean_mod
+! Renamed output variables
+! Add UTM coordinates to mask, 3d, xy, xz, yz output
+! 
+! 3337 2018-10-12 15:17:09Z kanani
 ! (from branch resler)
 ! Add biometeorology
@@ -330,11 +335,11 @@
     INTEGER(iwp), PARAMETER ::  dopr_norm_num = 7, dopts_num = 29, dots_max = 100
 
-    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
-         (/ 'wpt0  ', 'ws2   ', 'tsw2  ', 'ws3   ', 'ws2tsw', 'wstsw2',        &
-            'z_i   ' /)
-
-    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
-         (/ 'wpt0  ', 'w*2   ', 't*w2  ', 'w*3   ', 'w*2t*w', 'w*t*w2',        &
-            'z_i   ' /)
+    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
+         (/ 'wtheta0', 'ws2    ', 'tsw2   ', 'ws3    ', 'ws2tsw ', 'wstsw2 ',  &
+            'z_i    ' /)
+
+    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
+         (/ 'wtheta0', 'w*2    ', 't*w2   ', 'w*3    ', 'w*2t*w ', 'w*t*w2 ',  &
+            'z_i    ' /)
 
     CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
@@ -358,9 +363,9 @@
     CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
           (/ 'E            ', 'E*           ', 'dt           ',                &
-             'u*           ', 'th*          ', 'umax         ',                &
+             'us*          ', 'th*          ', 'umax         ',                &
              'vmax         ', 'wmax         ', 'div_new      ',                &
-             'div_old      ', 'z_i_wpt      ', 'z_i_pt       ',                &
-             'w*           ', 'w"pt"0       ', 'w"pt"        ',                &
-             'wpt          ', 'pt(0)        ', 'pt(z_mo)     ',                &
+             'div_old      ', 'zi_wtheta    ', 'zi_theta     ',                &
+             'w*           ', 'w"theta"0    ', 'w"theta"     ',                &
+             'wtheta       ', 'theta(0)     ', 'theta(z_mo)  ',                &
              'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
              'ol           ', 'q*           ', 'w"s"         ',                &
@@ -456,4 +461,10 @@
                     id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
 
+    INTEGER(iwp), DIMENSION(0:1,0:2) ::  id_var_eutm_3d, id_var_nutm_3d, &
+                                         id_var_eutm_xy, id_var_nutm_xy, &
+                                         id_var_eutm_xz, id_var_nutm_xz, &
+                                         id_var_eutm_yz, id_var_nutm_yz, &
+                                         id_var_eutm_ma, id_var_nutm_ma
+
     INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
     INTEGER ::  netcdf_deflate = 0      !< NetCDF compression, default: no
@@ -498,4 +509,7 @@
                    id_var_zu_mask, id_var_zw_mask, &
                    id_var_zusi_mask, id_var_zwwi_mask
+                   
+    INTEGER(iwp), DIMENSION(1:max_masks,0:1,0:2) ::  id_var_eutm_mask, &
+                                                     id_var_nutm_mask
 
     INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask
@@ -596,4 +610,7 @@
     USE land_surface_model_mod,                                                &
         ONLY: lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs
+
+    USE netcdf_data_input_mod,                                                 &
+        ONLY: init_model
 
     USE ocean_mod,                                                             &
@@ -692,4 +709,9 @@
 
     LOGICAL, SAVE ::  init_netcdf = .FALSE.                  !<
+
+    REAL(wp) ::  cos_ra                                      !< cosine of rotation_angle
+    REAL(wp) ::  shift_x                                     !< shift of x coordinate
+    REAL(wp) ::  shift_y                                     !< shift of y coordinate
+    REAL(wp) ::  sin_ra                                      !< sine of rotation_angle
 
     REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
@@ -901,4 +923,57 @@
                                   'meters', '', 486, 487, 000 )
 !
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_y_mask(mid,av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_xu_mask(mid,av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_y_mask(mid,av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_yv_mask(mid,av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
+!
 !--       In case of non-flat topography define 2d-arrays containing the height
 !--       information. Only for parallel netcdf output.
@@ -954,7 +1029,7 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',                &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       's', 'vpt'  )
+                       's', 'theta', 'thetal', 'thetav' )
 
                    grid_x = 'x'
@@ -1144,4 +1219,125 @@
 
 !
+!--       Write UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--          1D in case of no rotation
+             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--          x coordinates
+             ALLOCATE( netcdf_data(mask_size(mid,1)) )
+             DO  k = 0, 2
+!            
+!--             Scalar grid points
+                IF ( k == 0 )  THEN
+                   shift_x = 0.5
+!            
+!--             u grid points
+                ELSEIF ( k == 1 )  THEN
+                   shift_x = 0.0
+!            
+!--             v grid points
+                ELSEIF ( k == 2 )  THEN
+                   shift_x = 0.5
+                ENDIF
+
+                netcdf_data = init_model%origin_x + cos_ra                     &
+                       * ( mask_i_global(mid,:mask_size(mid,1)) + shift_x ) * dx
+
+                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
+                                        id_var_eutm_mask(mid,av,k), &
+                                        netcdf_data, start = (/ 1 /), &
+                                        count = (/ mask_size(mid,1) /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+             ENDDO
+             DEALLOCATE( netcdf_data )
+!
+!--          y coordinates
+             ALLOCATE( netcdf_data(mask_size(mid,2)) )
+             DO  k = 0, 2
+!
+!--             Scalar grid points
+                IF ( k == 0 )  THEN
+                   shift_y = 0.5
+!
+!--             u grid points
+                ELSEIF ( k == 1 )  THEN
+                   shift_y = 0.5
+!
+!--             v grid points
+                ELSEIF ( k == 2 )  THEN
+                   shift_y = 0.0
+                ENDIF
+
+                netcdf_data = init_model%origin_y + cos_ra                     &
+                       * ( mask_j_global(mid,:mask_size(mid,2)) + shift_y ) * dy
+
+                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
+                                        id_var_nutm_mask(mid,av,k), &
+                                        netcdf_data, start = (/ 1 /), &
+                                        count = (/ mask_size(mid,2) /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+             ENDDO
+             DEALLOCATE( netcdf_data )
+
+          ELSE
+!
+!--          2D in case of rotation
+             ALLOCATE( netcdf_data_2d(mask_size(mid,1),mask_size(mid,2)) )
+             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+
+             DO  k = 0, 2
+!            
+!--            Scalar grid points
+               IF ( k == 0 )  THEN
+                  shift_x = 0.5 ; shift_y = 0.5
+!            
+!--            u grid points
+               ELSEIF ( k == 1 )  THEN
+                  shift_x = 0.0 ; shift_y = 0.5
+!            
+!--            v grid points
+               ELSEIF ( k == 2 )  THEN
+                  shift_x = 0.5 ; shift_y = 0.0
+               ENDIF
+
+               DO  j = 1, mask_size(mid,2)
+                  DO  i = 1, mask_size(mid,1)
+                     netcdf_data_2d(i,j) = init_model%origin_x                 &
+                           + cos_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
+                           + sin_ra * ( mask_j_global(mid,j) + shift_y ) * dy
+                  ENDDO
+               ENDDO
+
+               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
+                                       id_var_eutm_mask(mid,av,k), &
+                                       netcdf_data_2d, start = (/ 1, 1 /), &
+                                       count = (/ mask_size(mid,1), &
+                                                  mask_size(mid,2) /) )
+               CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+               DO  j = 1, mask_size(mid,2)
+                  DO  i = 1, mask_size(mid,1)
+                     netcdf_data_2d(i,j) = init_model%origin_y                 &
+                           - sin_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
+                           + cos_ra * ( mask_j_global(mid,j) + shift_y ) * dy
+                  ENDDO
+               ENDDO
+             
+               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
+                                       id_var_nutm_mask(mid,av,k), &
+                                       netcdf_data_2d, start = (/ 1, 1 /), &
+                                       count = (/ mask_size(mid,1), &
+                                                  mask_size(mid,2) /) )
+               CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+             
+             ENDDO
+             DEALLOCATE( netcdf_data_2d )
+          ENDIF
+
+!
 !--       Write zu and zw data (vertical axes)
           ALLOCATE( netcdf_data(mask_size(mid,3)) )
@@ -1456,4 +1652,57 @@
                                   362, 363, 000 )
 !
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_y_3d(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_3d(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_xu_3d(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_3d(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_y_3d(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_3d(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_3d(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_yv_3d(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_3d(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av), id_dim_y_3d(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av), id_dim_y_3d(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_3d(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_3d(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_3d(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_3d(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_3d(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
+!
 !--       In case of non-flat topography define 2d-arrays containing the height
 !--       information. Only output 2d topography information in case of parallel
@@ -1506,7 +1755,7 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       'rho_ocean', 's', 'sa', 'vpt' )
+                       's', 'theta', 'thetal', 'thetav' )
 
                    grid_x = 'x'
@@ -1551,4 +1800,10 @@
                       CALL lsm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
                                                    grid_y, grid_z )
+                   ENDIF
+!
+!--                Check for ocean quantities
+                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                      CALL ocean_define_netcdf_grid( do3d(av,i), found,  &
+                                                     grid_x, grid_y, grid_z )
                    ENDIF
 
@@ -1729,4 +1984,122 @@
              DEALLOCATE( netcdf_data )
 
+!
+!--          Write UTM coordinates
+             IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--             1D in case of no rotation
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--             x coordinates
+                ALLOCATE( netcdf_data(0:nx) )
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5
+!               
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0
+!               
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5
+                   ENDIF
+                
+                   DO  i = 0, nx
+                     netcdf_data(i) = init_model%origin_x            &
+                                    + cos_ra * ( i + shift_x ) * dx
+                   ENDDO
+                
+                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ nx+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+!
+!--             y coordinates
+                ALLOCATE( netcdf_data(0:ny) )
+                DO  k = 0, 2
+!
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_y = 0.5
+!
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_y = 0.5
+!
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_y = 0.0
+                   ENDIF
+
+                   DO  j = 0, ny
+                      netcdf_data(j) = init_model%origin_y            &
+                                     + cos_ra * ( j + shift_y ) * dy
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+
+             ELSE
+!
+!--             2D in case of rotation
+                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+                
+                DO  k = 0, 2
+!               
+!--               Scalar grid points
+                  IF ( k == 0 )  THEN
+                     shift_x = 0.5 ; shift_y = 0.5
+!               
+!--               u grid points
+                  ELSEIF ( k == 1 )  THEN
+                     shift_x = 0.0 ; shift_y = 0.5
+!               
+!--               v grid points
+                  ELSEIF ( k == 2 )  THEN
+                     shift_x = 0.5 ; shift_y = 0.0
+                  ENDIF
+                
+                  DO  j = 0, ny
+                     DO  i = 0, nx
+                        netcdf_data_2d(i,j) = init_model%origin_x            &
+                                            + cos_ra * ( i + shift_x ) * dx  &
+                                            + sin_ra * ( j + shift_y ) * dy
+                     ENDDO
+                  ENDDO
+                
+                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(av,k),  &
+                                          netcdf_data_2d, start = (/ 1, 1 /),   &
+                                          count = (/ nx+1, ny+1 /) )
+                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+                
+                  DO  j = 0, ny
+                     DO  i = 0, nx
+                        netcdf_data_2d(i,j) = init_model%origin_y            &
+                                            - sin_ra * ( i + shift_x ) * dx  &
+                                            + cos_ra * ( j + shift_y ) * dy
+                     ENDDO
+                  ENDDO
+                
+                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(av,k),  &
+                                          netcdf_data_2d, start = (/ 1, 1 /),   &
+                                          count = (/ nx+1, ny+1 /) )
+                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+                
+                ENDDO
+                DEALLOCATE( netcdf_data_2d )
+             ENDIF
 !
 !--          Write zu and zw data (vertical axes)
@@ -2242,4 +2615,57 @@
                                   365, 366, 000 )
 !
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_y_xy(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xy(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_xu_xy(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xy(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_y_xy(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xy(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xy(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_yv_xy(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xy(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av), id_dim_y_xy(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av), id_dim_y_xy(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xy(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xy(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xy(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xy(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xy(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
+!
 !--       In case of non-flat topography define 2d-arrays containing the height
 !--       information. Only for parallel netcdf output.
@@ -2289,9 +2715,10 @@
 !
 !--                   Most variables are defined on the zu grid
-                      CASE ( 'e_xy', 'lpt_xy', 'nc_xy', 'nr_xy', 'p_xy',       &
-                             'pc_xy', 'pr_xy', 'prr_xy', 'pt_xy', 'q_xy',      &
+                      CASE ( 'e_xy', 'nc_xy', 'nr_xy', 'p_xy',                 &
+                             'pc_xy', 'pr_xy', 'prr_xy', 'q_xy',               &
                              'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
-                             'ql_vp_xy', 'qr_xy', 'qv_xy', 'rho_ocean_xy',     &
-                             's_xy', 'sa_xy', 'vpt_xy' )
+                             'ql_vp_xy', 'qr_xy', 'qv_xy',                     &
+                             's_xy',                                           &
+                             'theta_xy', 'thetal_xy', 'thetav_xy' )
 
                          grid_x = 'x'
@@ -2335,4 +2762,11 @@
                          ENDIF
 
+!
+!--                      Check for ocean quantities
+                         IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                            CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
+                                                           grid_x, grid_y,     &
+                                                           grid_z )
+                         ENDIF
 !
 !--                      Check for radiation quantities
@@ -2591,7 +3025,124 @@
 
              DEALLOCATE( netcdf_data )
-
-          ENDIF
-
+!
+!--          Write UTM coordinates
+             IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--             1D in case of no rotation
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--             x coordinates
+                ALLOCATE( netcdf_data(0:nx) )
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5
+!               
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0
+!               
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5
+                   ENDIF
+                
+                   DO  i = 0, nx
+                     netcdf_data(i) = init_model%origin_x            &
+                                    + cos_ra * ( i + shift_x ) * dx
+                   ENDDO
+                
+                   nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_eutm_xy(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ nx+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+!
+!--             y coordinates
+                ALLOCATE( netcdf_data(0:ny) )
+                DO  k = 0, 2
+!
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_y = 0.5
+!
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_y = 0.5
+!
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_y = 0.0
+                   ENDIF
+
+                   DO  j = 0, ny
+                      netcdf_data(j) = init_model%origin_y            &
+                                     + cos_ra * ( j + shift_y ) * dy
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_nutm_xy(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+
+             ELSE
+!
+!--             2D in case of rotation
+                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+                
+                DO  k = 0, 2
+!               
+!--               Scalar grid points
+                  IF ( k == 0 )  THEN
+                     shift_x = 0.5 ; shift_y = 0.5
+!               
+!--               u grid points
+                  ELSEIF ( k == 1 )  THEN
+                     shift_x = 0.0 ; shift_y = 0.5
+!               
+!--               v grid points
+                  ELSEIF ( k == 2 )  THEN
+                     shift_x = 0.5 ; shift_y = 0.0
+                  ENDIF
+                
+                  DO  j = 0, ny
+                     DO  i = 0, nx
+                        netcdf_data_2d(i,j) = init_model%origin_x            &
+                                            + cos_ra * ( i + shift_x ) * dx  &
+                                            + sin_ra * ( j + shift_y ) * dy
+                     ENDDO
+                  ENDDO
+                
+                  nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_eutm_xy(av,k),  &
+                                          netcdf_data_2d, start = (/ 1, 1 /),   &
+                                          count = (/ nx+1, ny+1 /) )
+                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+                
+                  DO  j = 0, ny
+                     DO  i = 0, nx
+                        netcdf_data_2d(i,j) = init_model%origin_y            &
+                                            - sin_ra * ( i + shift_x ) * dx  &
+                                            + cos_ra * ( j + shift_y ) * dy
+                     ENDDO
+                  ENDDO
+                
+                  nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_nutm_xy(av,k),  &
+                                          netcdf_data_2d, start = (/ 1, 1 /),   &
+                                          count = (/ nx+1, ny+1 /) )
+                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+                
+                ENDDO
+                DEALLOCATE( netcdf_data_2d )
+             ENDIF
+
+          ENDIF
 !
 !--       In case of non-flat topography write height information. Only for 
@@ -2985,4 +3536,57 @@
                                   NF90_DOUBLE, id_var_zw_xz(av), 'meters', '', &
                                   157, 158, 000 )
+!
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_y_xz(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_xu_xz(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_y_xz(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xz(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_yv_xz(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xz(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av), id_dim_y_xz(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av), id_dim_y_xz(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xz(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xz(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
 
           IF ( land_surface )  THEN
@@ -3011,8 +3615,9 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_xz', 'lpt_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz', &
-                          'pr_xz', 'prr_xz', 'pt_xz', 'q_xz', 'qc_xz',         &
+                   CASE ( 'e_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz',           &
+                          'pr_xz', 'prr_xz', 'q_xz', 'qc_xz',                  &
                           'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
-                          'qv_xz', 'rho_ocean_xz', 's_xz', 'sa_xz', 'vpt_xz' )
+                          'qv_xz', 's_xz',                                     &
+                          'theta_xz', 'thetal_xz', 'thetav_xz'                 )
 
                       grid_x = 'x'
@@ -3055,4 +3660,10 @@
                       ENDIF
 
+!
+!--                   Check for ocean quantities
+                      IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                         CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
+                                                        grid_x, grid_y, grid_z )
+                      ENDIF
 !
 !--                   Check for radiation quantities
@@ -3294,6 +3905,135 @@
              ENDIF
 
-
              DEALLOCATE( netcdf_data )
+!
+!--          Write UTM coordinates
+             IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--             1D in case of no rotation
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--             x coordinates
+                ALLOCATE( netcdf_data(0:nx) )
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5
+!               
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0
+!               
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5
+                   ENDIF
+                
+                   DO  i = 0, nx
+                     netcdf_data(i) = init_model%origin_x            &
+                                    + cos_ra * ( i + shift_x ) * dx
+                   ENDDO
+                
+                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_eutm_xz(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ nx+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+!
+!--             y coordinates
+                ALLOCATE( netcdf_data(1:ns) )
+                DO  k = 0, 2
+!
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_y = 0.5
+!
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_y = 0.5
+!
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_y = 0.0
+                   ENDIF
+
+                   DO  i = 1, ns
+                      IF( section(i,2) == -1 )  THEN
+                         netcdf_data(i) = -1.0_wp  ! section averaged along y
+                      ELSE
+                         netcdf_data(i) = init_model%origin_y &
+                                     + cos_ra * ( section(i,2) + shift_y ) * dy
+                      ENDIF
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_nutm_xz(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ns /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+
+             ELSE
+!
+!--             2D in case of rotation
+                ALLOCATE( netcdf_data_2d(0:nx,1:ns) )
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+                
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5 ; shift_y = 0.5
+!                  
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0 ; shift_y = 0.5
+!                  
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5 ; shift_y = 0.0
+                   ENDIF
+
+                   DO  j = 1, ns
+                      IF( section(j,2) == -1 )  THEN
+                         netcdf_data_2d(:,j) = -1.0_wp  ! section averaged along y
+                      ELSE
+                         DO  i = 0, nx
+                            netcdf_data_2d(i,j) = init_model%origin_x          &
+                                    + cos_ra * ( i + shift_x ) * dx            &
+                                    + sin_ra * ( section(j,2) + shift_y ) * dy
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                   
+                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_eutm_xz(av,k),  &
+                                           netcdf_data_2d, start = (/ 1, 1 /),   &
+                                           count = (/ nx+1, ns /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+                   
+                   DO  j = 1, ns
+                      IF( section(j,2) == -1 )  THEN
+                         netcdf_data_2d(:,j) = -1.0_wp  ! section averaged along y
+                      ELSE
+                         DO  i = 0, nx
+                            netcdf_data_2d(i,j) = init_model%origin_y          &
+                                    - sin_ra * ( i + shift_x ) * dx            &
+                                    + cos_ra * ( section(j,2) + shift_y ) * dy
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                   
+                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_nutm_xz(av,k),  &
+                                           netcdf_data_2d, start = (/ 1, 1 /),   &
+                                           count = (/ nx+1, ns /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+                
+                ENDDO
+                DEALLOCATE( netcdf_data_2d )
+             ENDIF
 
           ENDIF
@@ -3654,4 +4394,57 @@
                                   NF90_DOUBLE, id_var_zw_yz(av), 'meters', '', &
                                   196, 197, 000 )
+!
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_y_yz(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_yz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_xu_yz(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_yz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_y_yz(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_yz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_yz(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_yv_yz(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_yz(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av), id_dim_y_yz(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av), id_dim_y_yz(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_yz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_yz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_yz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_yz(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_yz(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
 
           IF ( land_surface )  THEN
@@ -3680,8 +4473,9 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_yz', 'lpt_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz', &
-                          'pr_yz','prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz', &
+                   CASE ( 'e_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz',           &
+                          'pr_yz','prr_yz', 'q_yz', 'qc_yz', 'ql_yz',          &
                           'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
-                          'rho_ocean_yz', 's_yz', 'sa_yz', 'vpt_yz' )
+                          's_yz',                                              &
+                          'theta_yz', 'thetal_yz', 'thetav_yz' )
 
                       grid_x = 'x'
@@ -3724,4 +4518,10 @@
                       ENDIF
 
+!
+!--                   Check for ocean quantities
+                      IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                         CALL ocean_define_netcdf_grid( do2d(av,i), found,     &
+                                                       grid_x, grid_y, grid_z )
+                      ENDIF
 !
 !--                   Check for radiation quantities
@@ -3951,4 +4751,134 @@
 
              DEALLOCATE( netcdf_data )
+!
+!--          Write UTM coordinates
+             IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--             1D in case of no rotation
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--             x coordinates
+                ALLOCATE( netcdf_data(1:ns) )
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5
+!               
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0
+!               
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5
+                   ENDIF
+                
+                   DO  i = 1, ns
+                      IF( section(i,3) == -1 )  THEN
+                         netcdf_data(i) = -1.0_wp  ! section averaged along x
+                      ELSE
+                         netcdf_data(i) = init_model%origin_x &
+                                     + cos_ra * ( section(i,3) + shift_x ) * dx
+                      ENDIF
+                   ENDDO
+                
+                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_eutm_yz(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ns /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+!
+!--             y coordinates
+                ALLOCATE( netcdf_data(0:ny) )
+                DO  k = 0, 2
+!
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_y = 0.5
+!
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_y = 0.5
+!
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_y = 0.0
+                   ENDIF
+
+                   DO  i = 0, ny
+                     netcdf_data(i) = init_model%origin_y            &
+                                    + cos_ra * ( i + shift_y ) * dy
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_nutm_yz(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+
+             ELSE
+!
+!--             2D in case of rotation
+                ALLOCATE( netcdf_data_2d(1:ns,0:ny) )
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5 ; shift_y = 0.5
+!                  
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0 ; shift_y = 0.5
+!                  
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5 ; shift_y = 0.0
+                   ENDIF
+
+                   DO  j = 0, ny
+                      DO  i = 1, ns
+                         IF( section(i,3) == -1 )  THEN
+                            netcdf_data_2d(i,:) = -1.0_wp !section averaged along x
+                         ELSE
+                            netcdf_data_2d(i,j) = init_model%origin_x          &
+                                    + cos_ra * ( section(i,3) + shift_x ) * dx &
+                                    + sin_ra * ( j + shift_y ) * dy
+                         ENDIF
+                      ENDDO
+                   ENDDO
+                   
+                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_eutm_yz(av,k),  &
+                                           netcdf_data_2d, start = (/ 1, 1 /),   &
+                                           count = (/ ns, ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+                   
+                   DO  j = 0, ny
+                      DO  i = 1, ns
+                         IF( section(i,3) == -1 )  THEN
+                            netcdf_data_2d(i,:) = -1.0_wp !section averaged along x
+                         ELSE
+                            netcdf_data_2d(i,j) = init_model%origin_y          &
+                                    - sin_ra * ( section(i,3) + shift_x ) * dx &
+                                    + cos_ra * ( j + shift_y ) * dy
+                         ENDIF
+                      ENDDO
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_nutm_yz(av,k),  &
+                                           netcdf_data_2d, start = (/ 1, 1 /),   &
+                                           count = (/ ns, ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+                
+                ENDDO
+                DEALLOCATE( netcdf_data_2d )
+             ENDIF
 
           ENDIF
@@ -4871,7 +5801,8 @@
 !
 !--             Most variables are defined on the zu levels
-                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       'rho_ocean', 's', 'sa', 'u', 'v', 'vpt' )
+                       'rho_sea_water', 's', 'sa', &
+                       'theta', 'thetal', 'thetav', 'u', 'v' )
 
                    grid_z = 'zu'
Index: /palm/trunk/SOURCE/ocean_mod.f90
===================================================================
--- /palm/trunk/SOURCE/ocean_mod.f90	(revision 3420)
+++ /palm/trunk/SOURCE/ocean_mod.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3381 2018-10-19 13:09:06Z raasch
 ! spin-up cooling for ocean surface implemented, see new parameter
 ! surface_cooling_spinup_time
@@ -674,5 +677,5 @@
     SELECT CASE ( TRIM( var ) )
 
-       CASE ( 'rho_ocean' )
+       CASE ( 'rho_sea_water' )
           unit = 'kg/m3'
 
@@ -723,5 +726,5 @@
           unit = dopr_unit
 
-       CASE ( 'rho_ocean' )
+       CASE ( 'rho_sea_water' )
           dopr_index(var_count) = 64
           dopr_unit             = 'kg/m3'
@@ -797,5 +800,6 @@
     SELECT CASE ( TRIM( var ) )
 
-       CASE ( 'rho_ocean', 'sa' )
+       CASE ( 'rho_sea_water', 'rho_sea_water_xy', 'rho_sea_water_xz', &
+              'rho_sea_water_yz', 'sa', 'sa_xy', 'sa_xz', 'sa_yz' )
           grid_x = 'x'
           grid_y = 'y'
@@ -846,5 +850,5 @@
        SELECT CASE ( TRIM( variable ) )
 
-          CASE ( 'rho_ocean' )
+          CASE ( 'rho_sea_water' )
              IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
                 ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -867,5 +871,5 @@
        SELECT CASE ( TRIM( variable ) )
 
-          CASE ( 'rho_ocean' )
+          CASE ( 'rho_sea_water' )
              IF ( ALLOCATED( rho_ocean_av ) )  THEN
                 DO  i = nxlg, nxrg
@@ -899,5 +903,5 @@
        SELECT CASE ( TRIM( variable ) )
 
-          CASE ( 'rho_ocean' )
+          CASE ( 'rho_sea_water' )
              IF ( ALLOCATED( rho_ocean_av ) ) THEN
                 DO  i = nxlg, nxrg
@@ -980,5 +984,5 @@
     SELECT CASE ( TRIM( variable ) )
 
-       CASE ( 'rho_ocean_xy', 'rho_ocean_xz', 'rho_ocean_yz' )
+       CASE ( 'rho_sea_water_xy', 'rho_sea_water_xz', 'rho_sea_water_yz' )
           IF ( av == 0 )  THEN
              to_be_resorted => rho_ocean
@@ -1073,5 +1077,5 @@
     SELECT CASE ( TRIM( variable ) )
 
-       CASE ( 'rho_ocean' )
+       CASE ( 'rho_sea_water' )
           IF ( av == 0 )  THEN
              to_be_resorted => rho_ocean
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 3420)
+++ /palm/trunk/SOURCE/parin.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added module for data output at surfaces
+! 
+! 3355 2018-10-16 14:03:34Z knoop
 ! - offline nesting separated from large-scale forcing module
 ! - top boundary condition for pressure in offline nesting changed
@@ -512,4 +515,7 @@
     USE statistics,                                                            &
         ONLY:  hom, hom_sum, pr_palm, region, statistic_regions
+
+    USE surface_output_mod,                                                    &
+        ONLY:  surface_output_parin
 
     USE synthetic_turbulence_generator_mod,                                    &
@@ -874,5 +880,11 @@
 !--       Check for module namelists and read them
           CALL lsm_parin
-          CALL bcm_parin
+!
+!--       Check if bulk cloud model is used and read namelist if required
+          CALL bcm_parin 
+!
+!--       Check if surface output is enabled and read 
+!--       &surface_output_parameters if required          
+          CALL surface_output_parin
           CALL usm_parin
           CALL spectra_parin
Index: /palm/trunk/SOURCE/print_1d.f90
===================================================================
--- /palm/trunk/SOURCE/print_1d.f90	(revision 3420)
+++ /palm/trunk/SOURCE/print_1d.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3241 2018-09-12 15:02:00Z raasch
 ! unused format statement removed
 ! 
@@ -181,9 +184,9 @@
 111 FORMAT (1X,131('-'))
 112 FORMAT (/)
-120 FORMAT ('   k     zu      u     du     v     dv     pt    dpt    ', &
+120 FORMAT ('   k     zu      u     du     v     dv     theta dtheta ', &
             ' e      Km    Kh     l      zu      k')
 121 FORMAT (1X,I4,1X,F7.1,1X,F6.2,1X,F5.2,1X,F6.2,1X,F5.2,2X,F6.2,1X,F5.2, &
             1X,F7.4,1X,F5.2,1X,F5.2,1X,F6.2,1X,F7.1,2X,I4)
-130 FORMAT ('   k     zw       w''pt''      wpt        w''u''       wu       ',&
+130 FORMAT ('   k     zw       w''theta''   wtheta     w''u''       wu       ',&
             '  w''v''       wv        zw      k')
 131 FORMAT (1X,I4,1X,F7.1,6(1X,E10.3),1X,F7.1,2X,I4)
Index: /palm/trunk/SOURCE/spectra_mod.f90
===================================================================
--- /palm/trunk/SOURCE/spectra_mod.f90	(revision 3420)
+++ /palm/trunk/SOURCE/spectra_mod.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3248 2018-09-14 09:42:06Z sward
 ! Minor formating changes
 ! 
@@ -622,5 +625,5 @@
           d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
        
-       CASE ( 'pt' )
+       CASE ( 'theta' )
           pr = 4
           d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
Index: /palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 3420)
+++ /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3337 2018-10-12 15:17:09Z kanani
 ! (from branch resler)
 ! Add biometeorology,
@@ -331,5 +334,5 @@
                 e_av = 0.0_wp
 
-             CASE ( 'lpt' )
+             CASE ( 'thetal' )
                 IF ( .NOT. ALLOCATED( lpt_av ) )  THEN
                    ALLOCATE( lpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -367,5 +370,5 @@
                 pr_av = 0.0_wp
 
-             CASE ( 'pt' )
+             CASE ( 'theta' )
                 IF ( .NOT. ALLOCATED( pt_av ) )  THEN
                    ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -457,5 +460,5 @@
                 u_av = 0.0_wp
 
-             CASE ( 'u*' )
+             CASE ( 'us*' )
                 IF ( .NOT. ALLOCATED( us_av ) )  THEN
                    ALLOCATE( us_av(nysg:nyng,nxlg:nxrg) )
@@ -469,5 +472,5 @@
                 v_av = 0.0_wp
 
-             CASE ( 'vpt' )
+             CASE ( 'thetav' )
                 IF ( .NOT. ALLOCATED( vpt_av ) )  THEN
                    ALLOCATE( vpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -609,5 +612,5 @@
              ENDIF
 
-          CASE ( 'lpt' )
+          CASE ( 'thetal' )
              IF ( ALLOCATED( lpt_av ) ) THEN
                 DO  i = nxlg, nxrg
@@ -712,5 +715,5 @@
              ENDIF
 
-          CASE ( 'pt' )
+          CASE ( 'theta' )
              IF ( ALLOCATED( pt_av ) ) THEN
                 IF ( .NOT. bulk_cloud_model ) THEN
@@ -1005,5 +1008,5 @@
              ENDIF
 
-          CASE ( 'u*' )
+          CASE ( 'us*' )
              IF ( ALLOCATED( us_av ) ) THEN   
                 DO  i = nxl, nxr
@@ -1044,5 +1047,5 @@
              ENDIF
 
-          CASE ( 'vpt' )
+          CASE ( 'thetav' )
              IF ( ALLOCATED( vpt_av ) ) THEN 
                 DO  i = nxlg, nxrg
Index: /palm/trunk/SOURCE/surface_output_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_output_mod.f90	(revision 3421)
+++ /palm/trunk/SOURCE/surface_output_mod.f90	(revision 3421)
@@ -0,0 +1,3286 @@
+!> @file surface_output_mod.f90
+!------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2018 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Add output variables
+! Write data into binary file
+! Initial implenentation from Klaus Ketelsen and Matthias Suehring
+!
+! 3420 2018-10-24 17:30:08Z gronemeier
+!
+! Description:
+! ------------
+!> Generate output for surface data.
+!
+!------------------------------------------------------------------------------!
+
+MODULE surface_output_mod
+
+   USE kinds
+
+   USE arrays_3d,                                                              &
+       ONLY:  zu, zw
+   USE control_parameters,                                                     &
+       ONLY:  surface_data_output
+   USE grid_variables,                                                         &
+       ONLY: dx,dy
+   USE indices,                                                                &
+       ONLY: nxl, nxr, nys, nyn, nzb, nzt
+   USE surface_mod,                                                            &
+       ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,                  &
+              surf_usm_h, surf_usm_v
+
+   IMPLICIT NONE
+
+   TYPE surf_out                !< data structure which contains all surfaces elements of all types on subdomain
+   
+      INTEGER(iwp) ::  ns       !< number of surface elements
+      INTEGER(iwp) ::  npoints  !< number of points  / vertices which define a surface element
+      
+      REAL(wp) ::  fillvalue = -9999.9_wp !< fillvalue for surface elements which are not defined
+      
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_out  !< output variables
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_av   !< variables used for averaging
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  points   !< points  / vertices of a surface element
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  polygons !< polygon data of a surface element
+   END TYPE surf_out
+   
+   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
+   
+   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
+
+   LOGICAL :: first_output(0:1) = .FALSE.                 !< true if first output was already called
+
+   REAL(wp) ::  averaging_interval_surf  = 9999999.9_wp   !< averaging interval
+   REAL(wp) ::  dt_dosurf = 9999999.9_wp                  !< time interval for instantaneous data output
+   REAL(wp) ::  dt_dosurf_av = 9999999.9_wp               !< time interval for averaged data output
+   REAL(wp) ::  skip_time_dosurf = 0.0_wp                 !< skip time for instantaneous data output
+   REAL(wp) ::  skip_time_dosurf_av = 0.0_wp              !< skip time for averaged data output
+   REAL(wp) ::  time_dosurf = 0.0_wp                      !< internal counter variable to check for instantaneous data output
+   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
+   
+   SAVE
+
+   PRIVATE
+
+   INTERFACE  surface_output
+      MODULE PROCEDURE surface_output
+   END INTERFACE  surface_output
+   
+   INTERFACE  surface_output_averaging
+      MODULE PROCEDURE surface_output_averaging
+   END INTERFACE  surface_output_averaging
+   
+   INTERFACE  surface_output_check_parameters
+      MODULE PROCEDURE surface_output_check_parameters
+   END INTERFACE  surface_output_check_parameters
+   
+   INTERFACE  surface_output_init
+      MODULE PROCEDURE surface_output_init
+   END INTERFACE  surface_output_init
+   
+   INTERFACE  surface_output_parin
+      MODULE PROCEDURE surface_output_parin
+   END INTERFACE  surface_output_parin
+
+!
+!--Public subroutines
+   PUBLIC surface_output, surface_output_averaging,                            &
+          surface_output_check_parameters, surface_output_init,                &
+          surface_output_parin
+!
+!--Public variables
+   PUBLIC average_count_surf, averaging_interval_surf, dt_dosurf, dt_dosurf_av,&
+          skip_time_dosurf, skip_time_dosurf_av, time_dosurf, time_dosurf_av
+
+ CONTAINS
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialization surface-data output data structure: calculation of vertices
+!> and polygon data for the surface elements, allocation of required arrays.
+!------------------------------------------------------------------------------!
+   SUBROUTINE surface_output_init
+      IMPLICIT NONE
+
+      INTEGER(iwp) ::  i     !< grid index in x-direction, also running variable for counting non-average data output
+      INTEGER(iwp) ::  ilen  !< string length
+      INTEGER(iwp) ::  j     !< grid index in y-direction, also running variable for counting average data output
+      INTEGER(iwp) ::  k     !< grid index in z-direction
+      INTEGER(iwp) ::  l     !< running index for surface-element orientation
+      INTEGER(iwp) ::  m     !< running index for surface elements
+      INTEGER(iwp) ::  n_out !< running index for number of output variables
+      
+     
+      INTEGER(iwp) ::  npg   !< counter variable for all surface elements ( or polygons )
+      
+      INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) ::  point_index !< dummy array used to check where the reference points for surface polygons are located
+           
+!
+!--   Determine the number of surface elements on subdomain
+      surfaces%ns = surf_def_h(0)%ns + surf_lsm_h%ns + surf_usm_h%ns           & !horizontal upward-facing
+                  + surf_def_h(1)%ns                                           & !horizontal downard-facing    
+                  + surf_def_v(0)%ns + surf_lsm_v(0)%ns + surf_usm_v(0)%ns     & !northward-facing
+                  + surf_def_v(1)%ns + surf_lsm_v(1)%ns + surf_usm_v(1)%ns     & !southward-facing   
+                  + surf_def_v(2)%ns + surf_lsm_v(2)%ns + surf_usm_v(2)%ns     & !westward-facing   
+                  + surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns       !eastward-facing   
+!
+!--   Allocate output variable and set to _FillValue attribute
+      ALLOCATE ( surfaces%var_out(1:surfaces%ns) )
+      surfaces%var_out = surfaces%fillvalue
+!
+!--   If there is an output of time average output variables, allocate the 
+!--   required array.
+      IF ( dosurf_no(1) > 0 )  THEN
+         ALLOCATE ( surfaces%var_av(1:surfaces%ns,1:dosurf_no(1)) )
+         surfaces%var_av = 0.0_wp
+      ENDIF
+!
+!--   Initialize points 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 ( point_index(k,j,i) < 0)  THEN
+               surfaces%npoints   = surfaces%npoints + 1    
+               point_index(k,j,i) = surfaces%npoints
+            ENDIF
+            IF ( point_index(k,j+1,i) < 0)  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k,j+1,i) = surfaces%npoints
+            ENDIF
+            IF ( point_index(k,j,i+1) < 0)  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k,j,i+1) = surfaces%npoints
+            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
+            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
+         ENDIF
+         IF ( point_index(k,j+1,i) < 0)  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j+1,i) = surfaces%npoints
+         ENDIF
+         IF ( point_index(k,j,i+1) < 0)  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j,i+1) = surfaces%npoints
+         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
+         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
+         ENDIF
+         IF ( point_index(k,j+1,i) < 0)  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j+1,i) = surfaces%npoints
+         ENDIF
+         IF ( point_index(k,j,i+1) < 0)  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j,i+1) = surfaces%npoints
+         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
+         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 == 2 )
+            j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 0 )
+            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
+            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
+            ENDIF
+!
+!--         Upper / lower right index for north- and south-facing surfaces
+            IF ( l == 0  .AND.  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
+               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
+               ENDIF
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .AND.  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
+               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
+               ENDIF
+            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 == 2 )
+            j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 0 )
+            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
+            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
+            ENDIF
+!
+!--         Upper / lower right index for north- and south-facing surfaces
+            IF ( l == 0  .AND.  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
+               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
+               ENDIF
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .AND.  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
+               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
+               ENDIF
+            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 == 2 )
+            j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 0 )
+            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
+            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
+            ENDIF
+!
+!--         Upper / lower right index for north- and south-facing surfaces
+            IF ( l == 0  .AND.  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
+               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
+               ENDIF
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .AND.  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
+               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
+               ENDIF
+            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)    )      
+!
+!--   After the number of vertices is counted, repeat the loops and define the
+!--   vertices. 
+!--   Start with the horizontal default surfaces
+      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) = surfaces%npoints
+               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+1,i) < 0)  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k,j+1,i) = surfaces%npoints
+               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
+            IF ( point_index(k,j,i+1) < 0)  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k,j,i+1) = surfaces%npoints
+               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) = surfaces%npoints
+               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
+            
+            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+1,i)
+            surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
+            surfaces%polygons(5,npg)   = point_index(k,j,i+1)  
+         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
+            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+1,i) < 0)  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j+1,i) = surfaces%npoints
+            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
+         IF ( point_index(k,j,i+1) < 0)  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j,i+1) = surfaces%npoints
+            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) = surfaces%npoints
+            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
+         
+         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+1,i)
+         surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
+         surfaces%polygons(5,npg)   = point_index(k,j,i+1)  
+      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
+            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+1,i) < 0)  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j+1,i) = surfaces%npoints
+            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
+         IF ( point_index(k,j,i+1) < 0)  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j,i+1) = surfaces%npoints
+            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) = surfaces%npoints
+            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
+         
+         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+1,i)
+         surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
+         surfaces%polygons(5,npg)   = point_index(k,j,i+1)    
+      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 == 2 )
+            j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 0 )
+            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
+               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 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
+               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  .AND.  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
+                  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+1,j,i+1) < 0)  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j,i+1) = surfaces%npoints
+                  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
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .AND.  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
+                  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 
+               IF ( point_index(k+1,j+1,i) < 0)  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j+1,i) = surfaces%npoints
+                  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
+            ENDIF          
+            
+            npg = npg + 1            
+            IF ( l == 0  .AND.  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 == 2 )
+            j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 0 )
+            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
+               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 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
+               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  .AND.  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
+                  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+1,j,i+1) < 0)  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j,i+1) = surfaces%npoints
+                  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
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .AND.  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
+                  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 
+               IF ( point_index(k+1,j+1,i) < 0)  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j+1,i) = surfaces%npoints
+                  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
+            ENDIF
+            
+            npg = npg + 1            
+            IF ( l == 0  .AND.  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 == 2 )
+            j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 0 )
+            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
+               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 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
+               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  .AND.  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
+                  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+1,j,i+1) < 0)  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j,i+1) = surfaces%npoints
+                  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
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .AND.  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
+                  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 
+               IF ( point_index(k+1,j+1,i) < 0)  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j+1,i) = surfaces%npoints
+                  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
+            ENDIF
+            
+            npg = npg + 1            
+            IF ( l == 0  .AND.  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
+!
+!--   Deallocate temporary dummy variable
+      DEALLOCATE ( point_index )
+!
+!--   Now, vertices as well as polygon data can be written to NetCDF files
+!--   ... 
+   END SUBROUTINE surface_output_init
+   
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Routine for controlling the data output. Surface data is collected from 
+!> different types of surfaces (default, natural, urban) and different 
+!> orientation and written to one 1D-output array. Further, NetCDF routines
+!> are called to write the surface data in the respective NetCDF files.
+!------------------------------------------------------------------------------!   
+   SUBROUTINE surface_output( av )
+   
+      USE control_parameters,                                                  &
+          ONLY:  io_blocks, io_group, simulated_time
+
+      USE pegrid,                                                              &
+          ONLY:  comm2d, ierr
+
+   
+      IMPLICIT NONE
+
+      CHARACTER(LEN=100) ::  trimvar = ' ' !< dummy for single output variable
+      
+      INTEGER(iwp) ::  av     !< id indicating average or non-average data output
+      INTEGER(iwp) ::  i      !< loop index
+      INTEGER(iwp) ::  l      !< running index for surface-element orientation
+      INTEGER(iwp) ::  m      !< running index for surface elements
+      INTEGER(iwp) ::  n_out  !< counter variables for surface output
+
+!
+!--   Return, if nothing to output
+      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%ns
+               WRITE ( 25+av )  surfaces%points
+               WRITE ( 25+av )  surfaces%polygons
+            ENDIF
+            CALL MPI_BARRIER( comm2d, ierr )
+            first_output(av) = .TRUE.
+         ENDDO
+      ENDIF
+!
+!--   Write time coordinate
+      DO  i = 0, io_blocks-1
+         IF ( i == io_group )  THEN
+            trimvar = 'time'
+            WRITE ( 25+av )  trimvar
+            WRITE ( 25+av )  simulated_time
+         ENDIF
+         CALL MPI_BARRIER( comm2d, ierr )
+      ENDDO
+
+
+      n_out = 0
+      DO  WHILE ( dosurf(av,n_out+1)(1:1) /= ' ' )
+
+         n_out   = n_out + 1
+         trimvar = TRIM( dosurf(av,n_out) )
+!
+!--      Set the output array to the _FillValue in case it is not 
+!--      defined for each type of surface.
+         surfaces%var_out = surfaces%fillvalue
+         SELECT CASE ( trimvar )
+
+            CASE ( 'us' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%us,               &
+                                             surf_def_h(1)%us,                 &
+                                             surf_lsm_h%us,                    &
+                                             surf_usm_h%us,                    &
+                                             surf_def_v(0)%us,                 &
+                                             surf_lsm_v(0)%us,                 &
+                                             surf_usm_v(0)%us,                 &
+                                             surf_def_v(1)%us,                 &
+                                             surf_lsm_v(1)%us,                 &
+                                             surf_usm_v(1)%us,                 &
+                                             surf_def_v(2)%us,                 &
+                                             surf_lsm_v(2)%us,                 &
+                                             surf_usm_v(2)%us,                 &
+                                             surf_def_v(3)%us,                 &
+                                             surf_lsm_v(3)%us,                 &
+                                             surf_usm_v(3)%us ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'ts' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%ts,               &
+                                               surf_def_h(1)%ts,               &
+                                               surf_lsm_h%ts,                  &
+                                               surf_usm_h%ts,                  &
+                                               surf_def_v(0)%ts,               &
+                                               surf_lsm_v(0)%ts,               &
+                                               surf_usm_v(0)%ts,               &
+                                               surf_def_v(1)%ts,               &
+                                               surf_lsm_v(1)%ts,               &
+                                               surf_usm_v(1)%ts,               &
+                                               surf_def_v(2)%ts,               &
+                                               surf_lsm_v(2)%ts,               &
+                                               surf_usm_v(2)%ts,               &
+                                               surf_def_v(3)%ts,               &
+                                               surf_lsm_v(3)%ts,               &
+                                               surf_usm_v(3)%ts ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'qs' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%qs,               &
+                                               surf_def_h(1)%qs,               &
+                                               surf_lsm_h%qs,                  &
+                                               surf_usm_h%qs,                  &
+                                               surf_def_v(0)%qs,               &
+                                               surf_lsm_v(0)%qs,               &
+                                               surf_usm_v(0)%qs,               &
+                                               surf_def_v(1)%qs,               &
+                                               surf_lsm_v(1)%qs,               &
+                                               surf_usm_v(1)%qs,               &
+                                               surf_def_v(2)%qs,               &
+                                               surf_lsm_v(2)%qs,               &
+                                               surf_usm_v(2)%qs,               &
+                                               surf_def_v(3)%qs,               &
+                                               surf_lsm_v(3)%qs,               &
+                                               surf_usm_v(3)%qs ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'ss' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%ss,               &
+                                               surf_def_h(1)%ss,               &
+                                               surf_lsm_h%ss,                  &
+                                               surf_usm_h%ss,                  &
+                                               surf_def_v(0)%ss,               &
+                                               surf_lsm_v(0)%ss,               &
+                                               surf_usm_v(0)%ss,               &
+                                               surf_def_v(1)%ss,               &
+                                               surf_lsm_v(1)%ss,               &
+                                               surf_usm_v(1)%ss,               &
+                                               surf_def_v(2)%ss,               &
+                                               surf_lsm_v(2)%ss,               &
+                                               surf_usm_v(2)%ss,               &
+                                               surf_def_v(3)%ss,               &
+                                               surf_lsm_v(3)%ss,               &
+                                               surf_usm_v(3)%ss ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'qcs' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%qcs,              &
+                                               surf_def_h(1)%qcs,              &
+                                               surf_lsm_h%qcs,                 &
+                                               surf_usm_h%qcs,                 &
+                                               surf_def_v(0)%qcs,              &
+                                               surf_lsm_v(0)%qcs,              &
+                                               surf_usm_v(0)%qcs,              &
+                                               surf_def_v(1)%qcs,              &
+                                               surf_lsm_v(1)%qcs,              &
+                                               surf_usm_v(1)%qcs,              &
+                                               surf_def_v(2)%qcs,              &
+                                               surf_lsm_v(2)%qcs,              &
+                                               surf_usm_v(2)%qcs,              &
+                                               surf_def_v(3)%qcs,              &
+                                               surf_lsm_v(3)%qcs,              &
+                                               surf_usm_v(3)%qcs ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'ncs' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%ncs,              &
+                                               surf_def_h(1)%ncs,              &
+                                               surf_lsm_h%ncs,                 &
+                                               surf_usm_h%ncs,                 &
+                                               surf_def_v(0)%ncs,              &
+                                               surf_lsm_v(0)%ncs,              &
+                                               surf_usm_v(0)%ncs,              &
+                                               surf_def_v(1)%ncs,              &
+                                               surf_lsm_v(1)%ncs,              &
+                                               surf_usm_v(1)%ncs,              &
+                                               surf_def_v(2)%ncs,              &
+                                               surf_lsm_v(2)%ncs,              &
+                                               surf_usm_v(2)%ncs,              &
+                                               surf_def_v(3)%ncs,              &
+                                               surf_lsm_v(3)%ncs,              &
+                                               surf_usm_v(3)%ncs ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'qrs' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%qrs,              &
+                                               surf_def_h(1)%qrs,              &
+                                               surf_lsm_h%qrs,                 &
+                                               surf_usm_h%qrs,                 &
+                                               surf_def_v(0)%qrs,              &
+                                               surf_lsm_v(0)%qrs,              &
+                                               surf_usm_v(0)%qrs,              &
+                                               surf_def_v(1)%qrs,              &
+                                               surf_lsm_v(1)%qrs,              &
+                                               surf_usm_v(1)%qrs,              &
+                                               surf_def_v(2)%qrs,              &
+                                               surf_lsm_v(2)%qrs,              &
+                                               surf_usm_v(2)%qrs,              &
+                                               surf_def_v(3)%qrs,              &
+                                               surf_lsm_v(3)%qrs,              &
+                                               surf_usm_v(3)%qrs ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'nrs' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%nrs,              &
+                                               surf_def_h(1)%nrs,              &
+                                               surf_lsm_h%nrs,                 &
+                                               surf_usm_h%nrs,                 &
+                                               surf_def_v(0)%nrs,              &
+                                               surf_lsm_v(0)%nrs,              &
+                                               surf_usm_v(0)%nrs,              &
+                                               surf_def_v(1)%nrs,              &
+                                               surf_lsm_v(1)%nrs,              &
+                                               surf_usm_v(1)%nrs,              &
+                                               surf_def_v(2)%nrs,              &
+                                               surf_lsm_v(2)%nrs,              &
+                                               surf_usm_v(2)%nrs,              &
+                                               surf_def_v(3)%nrs,              &
+                                               surf_lsm_v(3)%nrs,              &
+                                               surf_usm_v(3)%nrs ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'ol' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%ol,               &
+                                               surf_def_h(1)%ol,               &
+                                               surf_lsm_h%ol,                  &
+                                               surf_usm_h%ol,                  &
+                                               surf_def_v(0)%ol,               &
+                                               surf_lsm_v(0)%ol,               &
+                                               surf_usm_v(0)%ol,               &
+                                               surf_def_v(1)%ol,               &
+                                               surf_lsm_v(1)%ol,               &
+                                               surf_usm_v(1)%ol,               &
+                                               surf_def_v(2)%ol,               &
+                                               surf_lsm_v(2)%ol,               &
+                                               surf_usm_v(2)%ol,               &
+                                               surf_def_v(3)%ol,               &
+                                               surf_lsm_v(3)%ol,               &
+                                               surf_usm_v(3)%ol ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'z0' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%z0,               &
+                                               surf_def_h(1)%z0,               &
+                                               surf_lsm_h%z0,                  &
+                                               surf_usm_h%z0,                  &
+                                               surf_def_v(0)%z0,               &
+                                               surf_lsm_v(0)%z0,               &
+                                               surf_usm_v(0)%z0,               &
+                                               surf_def_v(1)%z0,               &
+                                               surf_lsm_v(1)%z0,               &
+                                               surf_usm_v(1)%z0,               &
+                                               surf_def_v(2)%z0,               &
+                                               surf_lsm_v(2)%z0,               &
+                                               surf_usm_v(2)%z0,               &
+                                               surf_def_v(3)%z0,               &
+                                               surf_lsm_v(3)%z0,               &
+                                               surf_usm_v(3)%z0 )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'z0h' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%z0h,              &
+                                               surf_def_h(1)%z0h,              &
+                                               surf_lsm_h%z0h,                 &
+                                               surf_usm_h%z0h,                 &
+                                               surf_def_v(0)%z0h,              &
+                                               surf_lsm_v(0)%z0h,              &
+                                               surf_usm_v(0)%z0h,              &
+                                               surf_def_v(1)%z0h,              &
+                                               surf_lsm_v(1)%z0h,              &
+                                               surf_usm_v(1)%z0h,              &
+                                               surf_def_v(2)%z0h,              &
+                                               surf_lsm_v(2)%z0h,              &
+                                               surf_usm_v(2)%z0h,              &
+                                               surf_def_v(3)%z0h,              &
+                                               surf_lsm_v(3)%z0h,              &
+                                               surf_usm_v(3)%z0h )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'z0q' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%z0q,              &
+                                               surf_def_h(1)%z0q,              &
+                                               surf_lsm_h%z0q,                 &
+                                               surf_usm_h%z0q,                 &
+                                               surf_def_v(0)%z0q,              &
+                                               surf_lsm_v(0)%z0q,              &
+                                               surf_usm_v(0)%z0q,              &
+                                               surf_def_v(1)%z0q,              &
+                                               surf_lsm_v(1)%z0q,              &
+                                               surf_usm_v(1)%z0q,              &
+                                               surf_def_v(2)%z0q,              &
+                                               surf_lsm_v(2)%z0q,              &
+                                               surf_usm_v(2)%z0q,              &
+                                               surf_def_v(3)%z0q,              &
+                                               surf_lsm_v(3)%z0q,              &
+                                               surf_usm_v(3)%z0q ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+
+               ENDIF
+
+            CASE ( 'theta1' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%pt1,              &
+                                               surf_def_h(1)%pt1,              &
+                                               surf_lsm_h%pt1,                 &
+                                               surf_usm_h%pt1,                 &
+                                               surf_def_v(0)%pt1,              &
+                                               surf_lsm_v(0)%pt1,              &
+                                               surf_usm_v(0)%pt1,              &
+                                               surf_def_v(1)%pt1,              &
+                                               surf_lsm_v(1)%pt1,              &
+                                               surf_usm_v(1)%pt1,              &
+                                               surf_def_v(2)%pt1,              &
+                                               surf_lsm_v(2)%pt1,              &
+                                               surf_usm_v(2)%pt1,              &
+                                               surf_def_v(3)%pt1,              &
+                                               surf_lsm_v(3)%pt1,              &
+                                               surf_usm_v(3)%pt1 )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'qv1' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%qv1,              &
+                                               surf_def_h(1)%qv1,              &
+                                               surf_lsm_h%qv1,                 &
+                                               surf_usm_h%qv1,                 &
+                                               surf_def_v(0)%qv1,              &
+                                               surf_lsm_v(0)%qv1,              &
+                                               surf_usm_v(0)%qv1,              &
+                                               surf_def_v(1)%qv1,              &
+                                               surf_lsm_v(1)%qv1,              &
+                                               surf_usm_v(1)%qv1,              &
+                                               surf_def_v(2)%qv1,              &
+                                               surf_lsm_v(2)%qv1,              &
+                                               surf_usm_v(2)%qv1,              &
+                                               surf_def_v(3)%qv1,              &
+                                               surf_lsm_v(3)%qv1,              &
+                                               surf_usm_v(3)%qv1 )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'thetav1' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%vpt1,             &
+                                               surf_def_h(1)%vpt1,             &
+                                               surf_lsm_h%vpt1,                &
+                                               surf_usm_h%vpt1,                &
+                                               surf_def_v(0)%vpt1,             &
+                                               surf_lsm_v(0)%vpt1,             &
+                                               surf_usm_v(0)%vpt1,             &
+                                               surf_def_v(1)%vpt1,             &
+                                               surf_lsm_v(1)%vpt1,             &
+                                               surf_usm_v(1)%vpt1,             &
+                                               surf_def_v(2)%vpt1,             &
+                                               surf_lsm_v(2)%vpt1,             &
+                                               surf_usm_v(2)%vpt1,             &
+                                               surf_def_v(3)%vpt1,             &
+                                               surf_lsm_v(3)%vpt1,             &
+                                               surf_usm_v(3)%vpt1 )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'usws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%usws,             &
+                                               surf_def_h(1)%usws,             &
+                                               surf_lsm_h%usws,                &
+                                               surf_usm_h%usws,                &
+                                               surf_def_v(0)%usws,             &
+                                               surf_lsm_v(0)%usws,             &
+                                               surf_usm_v(0)%usws,             &
+                                               surf_def_v(1)%usws,             &
+                                               surf_lsm_v(1)%usws,             &
+                                               surf_usm_v(1)%usws,             &
+                                               surf_def_v(2)%usws,             &
+                                               surf_lsm_v(2)%usws,             &
+                                               surf_usm_v(2)%usws,             &
+                                               surf_def_v(3)%usws,             &
+                                               surf_lsm_v(3)%usws,             &
+                                               surf_usm_v(3)%usws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'vsws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%vsws,             &
+                                               surf_def_h(1)%vsws,             &
+                                               surf_lsm_h%vsws,                &
+                                               surf_usm_h%vsws,                &
+                                               surf_def_v(0)%vsws,             &
+                                               surf_lsm_v(0)%vsws,             &
+                                               surf_usm_v(0)%vsws,             &
+                                               surf_def_v(1)%vsws,             &
+                                               surf_lsm_v(1)%vsws,             &
+                                               surf_usm_v(1)%vsws,             &
+                                               surf_def_v(2)%vsws,             &
+                                               surf_lsm_v(2)%vsws,             &
+                                               surf_usm_v(2)%vsws,             &
+                                               surf_def_v(3)%vsws,             &
+                                               surf_lsm_v(3)%vsws,             &
+                                               surf_usm_v(3)%vsws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'shf' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%shf,              &
+                                               surf_def_h(1)%shf,              &
+                                               surf_lsm_h%shf,                 &
+                                               surf_usm_h%shf,                 &
+                                               surf_def_v(0)%shf,              &
+                                               surf_lsm_v(0)%shf,              &
+                                               surf_usm_v(0)%shf,              &
+                                               surf_def_v(1)%shf,              &
+                                               surf_lsm_v(1)%shf,              &
+                                               surf_usm_v(1)%shf,              &
+                                               surf_def_v(2)%shf,              &
+                                               surf_lsm_v(2)%shf,              &
+                                               surf_usm_v(2)%shf,              &
+                                               surf_def_v(3)%shf,              &
+                                               surf_lsm_v(3)%shf,              &
+                                               surf_usm_v(3)%shf )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+               ENDIF
+
+            CASE ( 'qsws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%qsws,             &
+                                               surf_def_h(1)%qsws,             &
+                                               surf_lsm_h%qsws,                &
+                                               surf_usm_h%qsws,                &
+                                               surf_def_v(0)%qsws,             &
+                                               surf_lsm_v(0)%qsws,             &
+                                               surf_usm_v(0)%qsws,             &
+                                               surf_def_v(1)%qsws,             &
+                                               surf_lsm_v(1)%qsws,             &
+                                               surf_usm_v(1)%qsws,             &
+                                               surf_def_v(2)%qsws,             &
+                                               surf_lsm_v(2)%qsws,             &
+                                               surf_usm_v(2)%qsws,             &
+                                               surf_def_v(3)%qsws,             &
+                                               surf_lsm_v(3)%qsws,             &
+                                               surf_usm_v(3)%qsws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'ssws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%ssws,             &
+                                               surf_def_h(1)%ssws,             &
+                                               surf_lsm_h%ssws,                &
+                                               surf_usm_h%ssws,                &
+                                               surf_def_v(0)%ssws,             &
+                                               surf_lsm_v(0)%ssws,             &
+                                               surf_usm_v(0)%ssws,             &
+                                               surf_def_v(1)%ssws,             &
+                                               surf_lsm_v(1)%ssws,             &
+                                               surf_usm_v(1)%ssws,             &
+                                               surf_def_v(2)%ssws,             &
+                                               surf_lsm_v(2)%ssws,             &
+                                               surf_usm_v(2)%ssws,             &
+                                               surf_def_v(3)%ssws,             &
+                                               surf_lsm_v(3)%ssws,             &
+                                               surf_usm_v(3)%ssws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'qcsws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%qcsws,            &
+                                               surf_def_h(1)%qcsws,            &
+                                               surf_lsm_h%qcsws,               &
+                                               surf_usm_h%qcsws,               &
+                                               surf_def_v(0)%qcsws,            &
+                                               surf_lsm_v(0)%qcsws,            &
+                                               surf_usm_v(0)%qcsws,            &
+                                               surf_def_v(1)%qcsws,            &
+                                               surf_lsm_v(1)%qcsws,            &
+                                               surf_usm_v(1)%qcsws,            &
+                                               surf_def_v(2)%qcsws,            &
+                                               surf_lsm_v(2)%qcsws,            &
+                                               surf_usm_v(2)%qcsws,            &
+                                               surf_def_v(3)%qcsws,            &
+                                               surf_lsm_v(3)%qcsws,            &
+                                               surf_usm_v(3)%qcsws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'ncsws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%ncsws,            &
+                                               surf_def_h(1)%ncsws,            &
+                                               surf_lsm_h%ncsws,               &
+                                               surf_usm_h%ncsws,               &
+                                               surf_def_v(0)%ncsws,            &
+                                               surf_lsm_v(0)%ncsws,            &
+                                               surf_usm_v(0)%ncsws,            &
+                                               surf_def_v(1)%ncsws,            &
+                                               surf_lsm_v(1)%ncsws,            &
+                                               surf_usm_v(1)%ncsws,            &
+                                               surf_def_v(2)%ncsws,            &
+                                               surf_lsm_v(2)%ncsws,            &
+                                               surf_usm_v(2)%ncsws,            &
+                                               surf_def_v(3)%ncsws,            &
+                                               surf_lsm_v(3)%ncsws,            &
+                                               surf_usm_v(3)%ncsws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'qrsws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%qrsws,            &
+                                               surf_def_h(1)%qrsws,            &
+                                               surf_lsm_h%qrsws,               &
+                                               surf_usm_h%qrsws,               &
+                                               surf_def_v(0)%qrsws,            &
+                                               surf_lsm_v(0)%qrsws,            &
+                                               surf_usm_v(0)%qrsws,            &
+                                               surf_def_v(1)%qrsws,            &
+                                               surf_lsm_v(1)%qrsws,            &
+                                               surf_usm_v(1)%qrsws,            &
+                                               surf_def_v(2)%qrsws,            &
+                                               surf_lsm_v(2)%qrsws,            &
+                                               surf_usm_v(2)%qrsws,            &
+                                               surf_def_v(3)%qrsws,            &
+                                               surf_lsm_v(3)%qrsws,            &
+                                               surf_usm_v(3)%qrsws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'nrsws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%nrsws,            &
+                                               surf_def_h(1)%nrsws,            &
+                                               surf_lsm_h%nrsws,               &
+                                               surf_usm_h%nrsws,               &
+                                               surf_def_v(0)%nrsws,            &
+                                               surf_lsm_v(0)%nrsws,            &
+                                               surf_usm_v(0)%nrsws,            &
+                                               surf_def_v(1)%nrsws,            &
+                                               surf_lsm_v(1)%nrsws,            &
+                                               surf_usm_v(1)%nrsws,            &
+                                               surf_def_v(2)%nrsws,            &
+                                               surf_lsm_v(2)%nrsws,            &
+                                               surf_usm_v(2)%nrsws,            &
+                                               surf_def_v(3)%nrsws,            &
+                                               surf_lsm_v(3)%nrsws,            &
+                                               surf_usm_v(3)%nrsws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'sasws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%sasws,            &
+                                               surf_def_h(1)%sasws,            &
+                                               surf_lsm_h%sasws,               &
+                                               surf_usm_h%sasws,               &
+                                               surf_def_v(0)%sasws,            &
+                                               surf_lsm_v(0)%sasws,            &
+                                               surf_usm_v(0)%sasws,            &
+                                               surf_def_v(1)%sasws,            &
+                                               surf_lsm_v(1)%sasws,            &
+                                               surf_usm_v(1)%sasws,            &
+                                               surf_def_v(2)%sasws,            &
+                                               surf_lsm_v(2)%sasws,            &
+                                               surf_usm_v(2)%sasws,            &
+                                               surf_def_v(3)%sasws,            &
+                                               surf_lsm_v(3)%sasws,            &
+                                               surf_usm_v(3)%sasws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'q_surface' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%q_surface,        &
+                                               surf_def_h(1)%q_surface,        &
+                                               surf_lsm_h%q_surface,           &
+                                               surf_usm_h%q_surface,           &
+                                               surf_def_v(0)%q_surface,        &
+                                               surf_lsm_v(0)%q_surface,        &
+                                               surf_usm_v(0)%q_surface,        &
+                                               surf_def_v(1)%q_surface,        &
+                                               surf_lsm_v(1)%q_surface,        &
+                                               surf_usm_v(1)%q_surface,        &
+                                               surf_def_v(2)%q_surface,        &
+                                               surf_lsm_v(2)%q_surface,        &
+                                               surf_usm_v(2)%q_surface,        &
+                                               surf_def_v(3)%q_surface,        &
+                                               surf_lsm_v(3)%q_surface,        &
+                                               surf_usm_v(3)%q_surface )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'theta_surface' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%pt_surface,       &
+                                               surf_def_h(1)%pt_surface,       &
+                                               surf_lsm_h%pt_surface,          &
+                                               surf_usm_h%pt_surface,          &
+                                               surf_def_v(0)%pt_surface,       &
+                                               surf_lsm_v(0)%pt_surface,       &
+                                               surf_usm_v(0)%pt_surface,       &
+                                               surf_def_v(1)%pt_surface,       &
+                                               surf_lsm_v(1)%pt_surface,       &
+                                               surf_usm_v(1)%pt_surface,       &
+                                               surf_def_v(2)%pt_surface,       &
+                                               surf_lsm_v(2)%pt_surface,       &
+                                               surf_usm_v(2)%pt_surface,       &
+                                               surf_def_v(3)%pt_surface,       &
+                                               surf_lsm_v(3)%pt_surface,       &
+                                               surf_usm_v(3)%pt_surface )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'thetav_surface' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%vpt_surface,      &
+                                               surf_def_h(1)%vpt_surface,      &
+                                               surf_lsm_h%vpt_surface,         &
+                                               surf_usm_h%vpt_surface,         &
+                                               surf_def_v(0)%vpt_surface,      &
+                                               surf_lsm_v(0)%vpt_surface,      &
+                                               surf_usm_v(0)%vpt_surface,      &
+                                               surf_def_v(1)%vpt_surface,      &
+                                               surf_lsm_v(1)%vpt_surface,      &
+                                               surf_usm_v(1)%vpt_surface,      &
+                                               surf_def_v(2)%vpt_surface,      &
+                                               surf_lsm_v(2)%vpt_surface,      &
+                                               surf_usm_v(2)%vpt_surface,      &
+                                               surf_def_v(3)%vpt_surface,      &
+                                               surf_lsm_v(3)%vpt_surface,      &
+                                               surf_usm_v(3)%vpt_surface)
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'rad_net' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%rad_net,          &
+                                               surf_def_h(1)%rad_net,          &
+                                               surf_lsm_h%rad_net,             &
+                                               surf_usm_h%rad_net,             &
+                                               surf_def_v(0)%rad_net,          &
+                                               surf_lsm_v(0)%rad_net,          &
+                                               surf_usm_v(0)%rad_net,          &
+                                               surf_def_v(1)%rad_net,          &
+                                               surf_lsm_v(1)%rad_net,          &
+                                               surf_usm_v(1)%rad_net,          &
+                                               surf_def_v(2)%rad_net,          &
+                                               surf_lsm_v(2)%rad_net,          &
+                                               surf_usm_v(2)%rad_net,          &
+                                               surf_def_v(3)%rad_net,          &
+                                               surf_lsm_v(3)%rad_net,          &
+                                               surf_usm_v(3)%rad_net ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'rad_lw_in' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%rad_lw_in,        &
+                                               surf_def_h(1)%rad_lw_in,        &
+                                               surf_lsm_h%rad_lw_in,           &
+                                               surf_usm_h%rad_lw_in,           &
+                                               surf_def_v(0)%rad_lw_in,        &
+                                               surf_lsm_v(0)%rad_lw_in,        &
+                                               surf_usm_v(0)%rad_lw_in,        &
+                                               surf_def_v(1)%rad_lw_in,        &
+                                               surf_lsm_v(1)%rad_lw_in,        &
+                                               surf_usm_v(1)%rad_lw_in,        &
+                                               surf_def_v(2)%rad_lw_in,        &
+                                               surf_lsm_v(2)%rad_lw_in,        &
+                                               surf_usm_v(2)%rad_lw_in,        &
+                                               surf_def_v(3)%rad_lw_in,        &
+                                               surf_lsm_v(3)%rad_lw_in,        &
+                                               surf_usm_v(3)%rad_lw_in )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'rad_lw_out' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%rad_lw_out,       &
+                                               surf_def_h(1)%rad_lw_out,       &
+                                               surf_lsm_h%rad_lw_out,          &
+                                               surf_usm_h%rad_lw_out,          &
+                                               surf_def_v(0)%rad_lw_out,       &
+                                               surf_lsm_v(0)%rad_lw_out,       &
+                                               surf_usm_v(0)%rad_lw_out,       &
+                                               surf_def_v(1)%rad_lw_out,       &
+                                               surf_lsm_v(1)%rad_lw_out,       &
+                                               surf_usm_v(1)%rad_lw_out,       &
+                                               surf_def_v(2)%rad_lw_out,       &
+                                               surf_lsm_v(2)%rad_lw_out,       &
+                                               surf_usm_v(2)%rad_lw_out,       &
+                                               surf_def_v(3)%rad_lw_out,       &
+                                               surf_lsm_v(3)%rad_lw_out,       &
+                                               surf_usm_v(3)%rad_lw_out )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'rad_sw_in' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%rad_sw_in,        &
+                                               surf_def_h(1)%rad_sw_in,        &
+                                               surf_lsm_h%rad_sw_in,           &
+                                               surf_usm_h%rad_sw_in,           &
+                                               surf_def_v(0)%rad_sw_in,        &
+                                               surf_lsm_v(0)%rad_sw_in,        &
+                                               surf_usm_v(0)%rad_sw_in,        &
+                                               surf_def_v(1)%rad_sw_in,        &
+                                               surf_lsm_v(1)%rad_sw_in,        &
+                                               surf_usm_v(1)%rad_sw_in,        &
+                                               surf_def_v(2)%rad_sw_in,        &
+                                               surf_lsm_v(2)%rad_sw_in,        &
+                                               surf_usm_v(2)%rad_sw_in,        &
+                                               surf_def_v(3)%rad_sw_in,        &
+                                               surf_lsm_v(3)%rad_sw_in,        &
+                                               surf_usm_v(3)%rad_sw_in )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'rad_sw_out' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%rad_sw_out,       &
+                                               surf_def_h(1)%rad_sw_out,       &
+                                               surf_lsm_h%rad_sw_out,          &
+                                               surf_usm_h%rad_sw_out,          &
+                                               surf_def_v(0)%rad_sw_out,       &
+                                               surf_lsm_v(0)%rad_sw_out,       &
+                                               surf_usm_v(0)%rad_sw_out,       &
+                                               surf_def_v(1)%rad_sw_out,       &
+                                               surf_lsm_v(1)%rad_sw_out,       &
+                                               surf_usm_v(1)%rad_sw_out,       &
+                                               surf_def_v(2)%rad_sw_out,       &
+                                               surf_lsm_v(2)%rad_sw_out,       &
+                                               surf_usm_v(2)%rad_sw_out,       &
+                                               surf_def_v(3)%rad_sw_out,       &
+                                               surf_lsm_v(3)%rad_sw_out,       &
+                                               surf_usm_v(3)%rad_sw_out )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'ghf' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%ghf,              &
+                                               surf_def_h(1)%ghf,              &
+                                               surf_lsm_h%ghf,                 &
+                                               surf_usm_h%ghf,                 &
+                                               surf_def_v(0)%ghf,              &
+                                               surf_lsm_v(0)%ghf,              &
+                                               surf_usm_v(0)%ghf,              &
+                                               surf_def_v(1)%ghf,              &
+                                               surf_lsm_v(1)%ghf,              &
+                                               surf_usm_v(1)%ghf,              &
+                                               surf_def_v(2)%ghf,              &
+                                               surf_lsm_v(2)%ghf,              &
+                                               surf_usm_v(2)%ghf,              &
+                                               surf_def_v(3)%ghf,              &
+                                               surf_lsm_v(3)%ghf,              &
+                                               surf_usm_v(3)%ghf )      
+                                                                        ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+                      
+            CASE ( 'r_a' )                                                     
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%r_a,              &
+                                               surf_def_h(1)%r_a,              &
+                                               surf_lsm_h%r_a,                 &
+                                               surf_usm_h%r_a,                 &
+                                               surf_def_v(0)%r_a,              &
+                                               surf_lsm_v(0)%r_a,              &
+                                               surf_usm_v(0)%r_a,              &
+                                               surf_def_v(1)%r_a,              &
+                                               surf_lsm_v(1)%r_a,              &
+                                               surf_usm_v(1)%r_a,              &
+                                               surf_def_v(2)%r_a,              &
+                                               surf_lsm_v(2)%r_a,              &
+                                               surf_usm_v(2)%r_a,              &
+                                               surf_def_v(3)%r_a,              &
+                                               surf_lsm_v(3)%r_a,              &
+                                               surf_usm_v(3)%r_a )      
+                                                                      ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+                        
+            CASE ( 'r_soil' )                                                  
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%r_soil,           &
+                                               surf_def_h(1)%r_soil,           &
+                                               surf_lsm_h%r_soil,              &
+                                               surf_usm_h%r_soil,              &
+                                               surf_def_v(0)%r_soil,           &
+                                               surf_lsm_v(0)%r_soil,           &
+                                               surf_usm_v(0)%r_soil,           &
+                                               surf_def_v(1)%r_soil,           &
+                                               surf_lsm_v(1)%r_soil,           &
+                                               surf_usm_v(1)%r_soil,           &
+                                               surf_def_v(2)%r_soil,           &
+                                               surf_lsm_v(2)%r_soil,           &
+                                               surf_usm_v(2)%r_soil,           &
+                                               surf_def_v(3)%r_soil,           &
+                                               surf_lsm_v(3)%r_soil,           &
+                                               surf_usm_v(3)%r_soil ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'r_canopy' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%r_canopy,         &
+                                               surf_def_h(1)%r_canopy,         &
+                                               surf_lsm_h%r_canopy,            &
+                                               surf_usm_h%r_canopy,            &
+                                               surf_def_v(0)%r_canopy,         &
+                                               surf_lsm_v(0)%r_canopy,         &
+                                               surf_usm_v(0)%r_canopy,         &
+                                               surf_def_v(1)%r_canopy,         &
+                                               surf_lsm_v(1)%r_canopy,         &
+                                               surf_usm_v(1)%r_canopy,         &
+                                               surf_def_v(2)%r_canopy,         &
+                                               surf_lsm_v(2)%r_canopy,         &
+                                               surf_usm_v(2)%r_canopy,         &
+                                               surf_def_v(3)%r_canopy,         &
+                                               surf_lsm_v(3)%r_canopy,         &
+                                               surf_usm_v(3)%r_canopy ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+            CASE ( 'r_s' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_output_collect( surf_def_h(0)%r_s,              &
+                                               surf_def_h(1)%r_s,              &
+                                               surf_lsm_h%r_s,                 &
+                                               surf_usm_h%r_s,                 &
+                                               surf_def_v(0)%r_s,              &
+                                               surf_lsm_v(0)%r_s,              &
+                                               surf_usm_v(0)%r_s,              &
+                                               surf_def_v(1)%r_s,              &
+                                               surf_lsm_v(1)%r_s,              &
+                                               surf_usm_v(1)%r_s,              &
+                                               surf_def_v(2)%r_s,              &
+                                               surf_lsm_v(2)%r_s,              &
+                                               surf_usm_v(2)%r_s,              &
+                                               surf_def_v(3)%r_s,              &
+                                               surf_lsm_v(3)%r_s,              &
+                                               surf_usm_v(3)%r_s ) 
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+                                                        
+               ENDIF
+
+!
+!--            Add further variables:
+!--            'css', 'cssws', 'qsws_liq', 'qsws_soil', 'qsws_veg'
+
+         END SELECT
+!
+!--      Write to binary file: 
+!--      - surfaces%points ( 3, 1-npoints )
+!--      - surfaces%polygons ( 5, 1-ns )
+!--      - surfaces%var_out ( 1-ns, time )
+!--      - 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 )  trimvar
+               WRITE ( 25+av )  surfaces%var_out
+            ENDIF
+            CALL MPI_BARRIER( comm2d, ierr )
+         ENDDO
+
+      ENDDO
+      
+!
+!--      If averaged output was written to NetCDF file, set the counter to zero
+         IF ( av == 1 )  average_count_surf = 0
+             
+   END SUBROUTINE surface_output
+   
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Routine for controlling the data averaging. 
+!------------------------------------------------------------------------------!     
+   SUBROUTINE surface_output_averaging
+   
+      IMPLICIT NONE
+      
+      CHARACTER(LEN=100) ::  trimvar !< dummy variable for current output variable
+       
+      INTEGER(iwp) ::  l      !< running index for surface-element orientation
+      INTEGER(iwp) ::  m      !< running index for surface elements
+      INTEGER(iwp) ::  n_out  !< counter variables for surface output
+      
+      n_out = 0
+      DO  WHILE ( dosurf(1,n_out+1)(1:1) /= ' ' )
+         
+         n_out   = n_out + 1
+         trimvar = TRIM( dosurf(1,n_out) )
+
+         SELECT CASE ( trimvar )
+
+            CASE ( 'us' )
+               CALL surface_output_sum_up( surf_def_h(0)%us,                   &
+                                           surf_def_h(1)%us,                   &
+                                           surf_lsm_h%us,                      &
+                                           surf_usm_h%us,                      &
+                                           surf_def_v(0)%us,                   &
+                                           surf_lsm_v(0)%us,                   &
+                                           surf_usm_v(0)%us,                   &
+                                           surf_def_v(1)%us,                   &
+                                           surf_lsm_v(1)%us,                   &
+                                           surf_usm_v(1)%us,                   &
+                                           surf_def_v(2)%us,                   &
+                                           surf_lsm_v(2)%us,                   &
+                                           surf_usm_v(2)%us,                   &
+                                           surf_def_v(3)%us,                   &
+                                           surf_lsm_v(3)%us,                   &
+                                           surf_usm_v(3)%us, n_out ) 
+
+            CASE ( 'ts' )
+               CALL surface_output_sum_up( surf_def_h(0)%ts,                   &
+                                           surf_def_h(1)%ts,                   &
+                                           surf_lsm_h%ts,                      &
+                                           surf_usm_h%ts,                      &
+                                           surf_def_v(0)%ts,                   &
+                                           surf_lsm_v(0)%ts,                   &
+                                           surf_usm_v(0)%ts,                   &
+                                           surf_def_v(1)%ts,                   &
+                                           surf_lsm_v(1)%ts,                   &
+                                           surf_usm_v(1)%ts,                   &
+                                           surf_def_v(2)%ts,                   &
+                                           surf_lsm_v(2)%ts,                   &
+                                           surf_usm_v(2)%ts,                   &
+                                           surf_def_v(3)%ts,                   &
+                                           surf_lsm_v(3)%ts,                   &
+                                           surf_usm_v(3)%ts, n_out ) 
+
+            CASE ( 'qs' )
+               CALL surface_output_sum_up( surf_def_h(0)%qs,                   &
+                                           surf_def_h(1)%qs,                   &
+                                           surf_lsm_h%qs,                      &
+                                           surf_usm_h%qs,                      &
+                                           surf_def_v(0)%qs,                   &
+                                           surf_lsm_v(0)%qs,                   &
+                                           surf_usm_v(0)%qs,                   &
+                                           surf_def_v(1)%qs,                   &
+                                           surf_lsm_v(1)%qs,                   &
+                                           surf_usm_v(1)%qs,                   &
+                                           surf_def_v(2)%qs,                   &
+                                           surf_lsm_v(2)%qs,                   &
+                                           surf_usm_v(2)%qs,                   &
+                                           surf_def_v(3)%qs,                   &
+                                           surf_lsm_v(3)%qs,                   &
+                                           surf_usm_v(3)%qs, n_out ) 
+
+            CASE ( 'ss' )
+               CALL surface_output_sum_up( surf_def_h(0)%ss,                   &
+                                           surf_def_h(1)%ss,                   &
+                                           surf_lsm_h%ss,                      &
+                                           surf_usm_h%ss,                      &
+                                           surf_def_v(0)%ss,                   &
+                                           surf_lsm_v(0)%ss,                   &
+                                           surf_usm_v(0)%ss,                   &
+                                           surf_def_v(1)%ss,                   &
+                                           surf_lsm_v(1)%ss,                   &
+                                           surf_usm_v(1)%ss,                   &
+                                           surf_def_v(2)%ss,                   &
+                                           surf_lsm_v(2)%ss,                   &
+                                           surf_usm_v(2)%ss,                   &
+                                           surf_def_v(3)%ss,                   &
+                                           surf_lsm_v(3)%ss,                   &
+                                           surf_usm_v(3)%ss, n_out ) 
+
+            CASE ( 'qcs' )
+               CALL surface_output_sum_up( surf_def_h(0)%qcs,                  &
+                                           surf_def_h(1)%qcs,                  &
+                                           surf_lsm_h%qcs,                     &
+                                           surf_usm_h%qcs,                     &
+                                           surf_def_v(0)%qcs,                  &
+                                           surf_lsm_v(0)%qcs,                  &
+                                           surf_usm_v(0)%qcs,                  &
+                                           surf_def_v(1)%qcs,                  &
+                                           surf_lsm_v(1)%qcs,                  &
+                                           surf_usm_v(1)%qcs,                  &
+                                           surf_def_v(2)%qcs,                  &
+                                           surf_lsm_v(2)%qcs,                  &
+                                           surf_usm_v(2)%qcs,                  &
+                                           surf_def_v(3)%qcs,                  &
+                                           surf_lsm_v(3)%qcs,                  &
+                                           surf_usm_v(3)%qcs, n_out ) 
+
+            CASE ( 'ncs' )
+               CALL surface_output_sum_up( surf_def_h(0)%ncs,                  &
+                                           surf_def_h(1)%ncs,                  &
+                                           surf_lsm_h%ncs,                     &
+                                           surf_usm_h%ncs,                     &
+                                           surf_def_v(0)%ncs,                  &
+                                           surf_lsm_v(0)%ncs,                  &
+                                           surf_usm_v(0)%ncs,                  &
+                                           surf_def_v(1)%ncs,                  &
+                                           surf_lsm_v(1)%ncs,                  &
+                                           surf_usm_v(1)%ncs,                  &
+                                           surf_def_v(2)%ncs,                  &
+                                           surf_lsm_v(2)%ncs,                  &
+                                           surf_usm_v(2)%ncs,                  &
+                                           surf_def_v(3)%ncs,                  &
+                                           surf_lsm_v(3)%ncs,                  &
+                                           surf_usm_v(3)%ncs, n_out ) 
+
+            CASE ( 'qrs' )
+               CALL surface_output_sum_up( surf_def_h(0)%qrs,                  &
+                                           surf_def_h(1)%qrs,                  &
+                                           surf_lsm_h%qrs,                     &
+                                           surf_usm_h%qrs,                     &
+                                           surf_def_v(0)%qrs,                  &
+                                           surf_lsm_v(0)%qrs,                  &
+                                           surf_usm_v(0)%qrs,                  &
+                                           surf_def_v(1)%qrs,                  &
+                                           surf_lsm_v(1)%qrs,                  &
+                                           surf_usm_v(1)%qrs,                  &
+                                           surf_def_v(2)%qrs,                  &
+                                           surf_lsm_v(2)%qrs,                  &
+                                           surf_usm_v(2)%qrs,                  &
+                                           surf_def_v(3)%qrs,                  &
+                                           surf_lsm_v(3)%qrs,                  &
+                                           surf_usm_v(3)%qrs, n_out ) 
+
+            CASE ( 'nrs' )
+               CALL surface_output_sum_up( surf_def_h(0)%nrs,                  &
+                                           surf_def_h(1)%nrs,                  &
+                                           surf_lsm_h%nrs,                     &
+                                           surf_usm_h%nrs,                     &
+                                           surf_def_v(0)%nrs,                  &
+                                           surf_lsm_v(0)%nrs,                  &
+                                           surf_usm_v(0)%nrs,                  &
+                                           surf_def_v(1)%nrs,                  &
+                                           surf_lsm_v(1)%nrs,                  &
+                                           surf_usm_v(1)%nrs,                  &
+                                           surf_def_v(2)%nrs,                  &
+                                           surf_lsm_v(2)%nrs,                  &
+                                           surf_usm_v(2)%nrs,                  &
+                                           surf_def_v(3)%nrs,                  &
+                                           surf_lsm_v(3)%nrs,                  &
+                                           surf_usm_v(3)%nrs, n_out ) 
+
+            CASE ( 'ol' )
+               CALL surface_output_sum_up( surf_def_h(0)%ol,                   &
+                                           surf_def_h(1)%ol,                   &
+                                           surf_lsm_h%ol,                      &
+                                           surf_usm_h%ol,                      &
+                                           surf_def_v(0)%ol,                   &
+                                           surf_lsm_v(0)%ol,                   &
+                                           surf_usm_v(0)%ol,                   &
+                                           surf_def_v(1)%ol,                   &
+                                           surf_lsm_v(1)%ol,                   &
+                                           surf_usm_v(1)%ol,                   &
+                                           surf_def_v(2)%ol,                   &
+                                           surf_lsm_v(2)%ol,                   &
+                                           surf_usm_v(2)%ol,                   &
+                                           surf_def_v(3)%ol,                   &
+                                           surf_lsm_v(3)%ol,                   &
+                                           surf_usm_v(3)%ol, n_out ) 
+
+            CASE ( 'z0' )
+               CALL surface_output_sum_up( surf_def_h(0)%z0,                   &
+                                           surf_def_h(1)%z0,                   &
+                                           surf_lsm_h%z0,                      &
+                                           surf_usm_h%z0,                      &
+                                           surf_def_v(0)%z0,                   &
+                                           surf_lsm_v(0)%z0,                   &
+                                           surf_usm_v(0)%z0,                   &
+                                           surf_def_v(1)%z0,                   &
+                                           surf_lsm_v(1)%z0,                   &
+                                           surf_usm_v(1)%z0,                   &
+                                           surf_def_v(2)%z0,                   &
+                                           surf_lsm_v(2)%z0,                   &
+                                           surf_usm_v(2)%z0,                   &
+                                           surf_def_v(3)%z0,                   &
+                                           surf_lsm_v(3)%z0,                   &
+                                           surf_usm_v(3)%z0, n_out ) 
+
+            CASE ( 'z0h' )
+               CALL surface_output_sum_up( surf_def_h(0)%z0h,                  &
+                                           surf_def_h(1)%z0h,                  &
+                                           surf_lsm_h%z0h,                     &
+                                           surf_usm_h%z0h,                     &
+                                           surf_def_v(0)%z0h,                  &
+                                           surf_lsm_v(0)%z0h,                  &
+                                           surf_usm_v(0)%z0h,                  &
+                                           surf_def_v(1)%z0h,                  &
+                                           surf_lsm_v(1)%z0h,                  &
+                                           surf_usm_v(1)%z0h,                  &
+                                           surf_def_v(2)%z0h,                  &
+                                           surf_lsm_v(2)%z0h,                  &
+                                           surf_usm_v(2)%z0h,                  &
+                                           surf_def_v(3)%z0h,                  &
+                                           surf_lsm_v(3)%z0h,                  &
+                                           surf_usm_v(3)%z0h, n_out )          
+
+            CASE ( 'z0q' )
+               CALL surface_output_sum_up( surf_def_h(0)%z0q,                  &
+                                           surf_def_h(1)%z0q,                  &
+                                           surf_lsm_h%z0q,                     &
+                                           surf_usm_h%z0q,                     &
+                                           surf_def_v(0)%z0q,                  &
+                                           surf_lsm_v(0)%z0q,                  &
+                                           surf_usm_v(0)%z0q,                  &
+                                           surf_def_v(1)%z0q,                  &
+                                           surf_lsm_v(1)%z0q,                  &
+                                           surf_usm_v(1)%z0q,                  &
+                                           surf_def_v(2)%z0q,                  &
+                                           surf_lsm_v(2)%z0q,                  &
+                                           surf_usm_v(2)%z0q,                  &
+                                           surf_def_v(3)%z0q,                  &
+                                           surf_lsm_v(3)%z0q,                  &
+                                           surf_usm_v(3)%z0q, n_out ) 
+
+            CASE ( 'theta1' )                                                  
+               CALL surface_output_sum_up( surf_def_h(0)%pt1,                  &
+                                           surf_def_h(1)%pt1,                  &
+                                           surf_lsm_h%pt1,                     &
+                                           surf_usm_h%pt1,                     &
+                                           surf_def_v(0)%pt1,                  &
+                                           surf_lsm_v(0)%pt1,                  &
+                                           surf_usm_v(0)%pt1,                  &
+                                           surf_def_v(1)%pt1,                  &
+                                           surf_lsm_v(1)%pt1,                  &
+                                           surf_usm_v(1)%pt1,                  &
+                                           surf_def_v(2)%pt1,                  &
+                                           surf_lsm_v(2)%pt1,                  &
+                                           surf_usm_v(2)%pt1,                  &
+                                           surf_def_v(3)%pt1,                  &
+                                           surf_lsm_v(3)%pt1,                  &
+                                           surf_usm_v(3)%pt1, n_out )          
+                                                                               
+            CASE ( 'qv1' )                                                     
+               CALL surface_output_sum_up( surf_def_h(0)%qv1,                  &
+                                           surf_def_h(1)%qv1,                  &
+                                           surf_lsm_h%qv1,                     &
+                                           surf_usm_h%qv1,                     &
+                                           surf_def_v(0)%qv1,                  &
+                                           surf_lsm_v(0)%qv1,                  &
+                                           surf_usm_v(0)%qv1,                  &
+                                           surf_def_v(1)%qv1,                  &
+                                           surf_lsm_v(1)%qv1,                  &
+                                           surf_usm_v(1)%qv1,                  &
+                                           surf_def_v(2)%qv1,                  &
+                                           surf_lsm_v(2)%qv1,                  &
+                                           surf_usm_v(2)%qv1,                  &
+                                           surf_def_v(3)%qv1,                  &
+                                           surf_lsm_v(3)%qv1,                  &
+                                           surf_usm_v(3)%qv1, n_out ) 
+
+            CASE ( 'thetav1' )
+               CALL surface_output_sum_up( surf_def_h(0)%vpt1,                 &
+                                           surf_def_h(1)%vpt1,                 &
+                                           surf_lsm_h%vpt1,                    &
+                                           surf_usm_h%vpt1,                    &
+                                           surf_def_v(0)%vpt1,                 &
+                                           surf_lsm_v(0)%vpt1,                 &
+                                           surf_usm_v(0)%vpt1,                 &
+                                           surf_def_v(1)%vpt1,                 &
+                                           surf_lsm_v(1)%vpt1,                 &
+                                           surf_usm_v(1)%vpt1,                 &
+                                           surf_def_v(2)%vpt1,                 &
+                                           surf_lsm_v(2)%vpt1,                 &
+                                           surf_usm_v(2)%vpt1,                 &
+                                           surf_def_v(3)%vpt1,                 &
+                                           surf_lsm_v(3)%vpt1,                 &
+                                           surf_usm_v(3)%vpt1, n_out ) 
+
+            CASE ( 'usws' )
+               CALL surface_output_sum_up( surf_def_h(0)%usws,                 &
+                                           surf_def_h(1)%usws,                 &
+                                           surf_lsm_h%usws,                    &
+                                           surf_usm_h%usws,                    &
+                                           surf_def_v(0)%usws,                 &
+                                           surf_lsm_v(0)%usws,                 &
+                                           surf_usm_v(0)%usws,                 &
+                                           surf_def_v(1)%usws,                 &
+                                           surf_lsm_v(1)%usws,                 &
+                                           surf_usm_v(1)%usws,                 &
+                                           surf_def_v(2)%usws,                 &
+                                           surf_lsm_v(2)%usws,                 &
+                                           surf_usm_v(2)%usws,                 &
+                                           surf_def_v(3)%usws,                 &
+                                           surf_lsm_v(3)%usws,                 &
+                                           surf_usm_v(3)%usws, n_out ) 
+
+            CASE ( 'vsws' )
+               CALL surface_output_sum_up( surf_def_h(0)%vsws,                 &
+                                           surf_def_h(1)%vsws,                 &
+                                           surf_lsm_h%vsws,                    &
+                                           surf_usm_h%vsws,                    &
+                                           surf_def_v(0)%vsws,                 &
+                                           surf_lsm_v(0)%vsws,                 &
+                                           surf_usm_v(0)%vsws,                 &
+                                           surf_def_v(1)%vsws,                 &
+                                           surf_lsm_v(1)%vsws,                 &
+                                           surf_usm_v(1)%vsws,                 &
+                                           surf_def_v(2)%vsws,                 &
+                                           surf_lsm_v(2)%vsws,                 &
+                                           surf_usm_v(2)%vsws,                 &
+                                           surf_def_v(3)%vsws,                 &
+                                           surf_lsm_v(3)%vsws,                 &
+                                           surf_usm_v(3)%vsws, n_out ) 
+
+            CASE ( 'shf' )
+               CALL surface_output_sum_up( surf_def_h(0)%shf,                  &
+                                           surf_def_h(1)%shf,                  &
+                                           surf_lsm_h%shf,                     &
+                                           surf_usm_h%shf,                     &  
+                                           surf_def_v(0)%shf,                  &
+                                           surf_lsm_v(0)%shf,                  &
+                                           surf_usm_v(0)%shf,                  &
+                                           surf_def_v(1)%shf,                  &
+                                           surf_lsm_v(1)%shf,                  &
+                                           surf_usm_v(1)%shf,                  &
+                                           surf_def_v(2)%shf,                  &
+                                           surf_lsm_v(2)%shf,                  &
+                                           surf_usm_v(2)%shf,                  &
+                                           surf_def_v(3)%shf,                  &
+                                           surf_lsm_v(3)%shf,                  &
+                                           surf_usm_v(3)%shf, n_out )
+
+            CASE ( 'qsws' )
+               CALL surface_output_sum_up( surf_def_h(0)%qsws,                 &
+                                           surf_def_h(1)%qsws,                 &
+                                           surf_lsm_h%qsws,                    &
+                                           surf_usm_h%qsws,                    &
+                                           surf_def_v(0)%qsws,                 &
+                                           surf_lsm_v(0)%qsws,                 &
+                                           surf_usm_v(0)%qsws,                 &
+                                           surf_def_v(1)%qsws,                 &
+                                           surf_lsm_v(1)%qsws,                 &
+                                           surf_usm_v(1)%qsws,                 &
+                                           surf_def_v(2)%qsws,                 &
+                                           surf_lsm_v(2)%qsws,                 &
+                                           surf_usm_v(2)%qsws,                 &
+                                           surf_def_v(3)%qsws,                 &
+                                           surf_lsm_v(3)%qsws,                 &
+                                           surf_usm_v(3)%qsws, n_out ) 
+
+            CASE ( 'ssws' )
+               CALL surface_output_sum_up( surf_def_h(0)%ssws,                 &
+                                           surf_def_h(1)%ssws,                 &
+                                           surf_lsm_h%ssws,                    &
+                                           surf_usm_h%ssws,                    &
+                                           surf_def_v(0)%ssws,                 &
+                                           surf_lsm_v(0)%ssws,                 &
+                                           surf_usm_v(0)%ssws,                 &
+                                           surf_def_v(1)%ssws,                 &
+                                           surf_lsm_v(1)%ssws,                 &
+                                           surf_usm_v(1)%ssws,                 &
+                                           surf_def_v(2)%ssws,                 &
+                                           surf_lsm_v(2)%ssws,                 &
+                                           surf_usm_v(2)%ssws,                 &
+                                           surf_def_v(3)%ssws,                 &
+                                           surf_lsm_v(3)%ssws,                 &
+                                           surf_usm_v(3)%ssws, n_out ) 
+
+            CASE ( 'qcsws' )
+               CALL surface_output_sum_up( surf_def_h(0)%qcsws,                &
+                                           surf_def_h(1)%qcsws,                &
+                                           surf_lsm_h%qcsws,                   &
+                                           surf_usm_h%qcsws,                   &
+                                           surf_def_v(0)%qcsws,                &
+                                           surf_lsm_v(0)%qcsws,                &
+                                           surf_usm_v(0)%qcsws,                &
+                                           surf_def_v(1)%qcsws,                &
+                                           surf_lsm_v(1)%qcsws,                &
+                                           surf_usm_v(1)%qcsws,                &
+                                           surf_def_v(2)%qcsws,                &
+                                           surf_lsm_v(2)%qcsws,                &
+                                           surf_usm_v(2)%qcsws,                &
+                                           surf_def_v(3)%qcsws,                &
+                                           surf_lsm_v(3)%qcsws,                &
+                                           surf_usm_v(3)%qcsws, n_out ) 
+
+            CASE ( 'ncsws' )
+               CALL surface_output_sum_up( surf_def_h(0)%ncsws,                &
+                                           surf_def_h(1)%ncsws,                &
+                                           surf_lsm_h%ncsws,                   &
+                                           surf_usm_h%ncsws,                   &
+                                           surf_def_v(0)%ncsws,                &
+                                           surf_lsm_v(0)%ncsws,                &
+                                           surf_usm_v(0)%ncsws,                &
+                                           surf_def_v(1)%ncsws,                &
+                                           surf_lsm_v(1)%ncsws,                &
+                                           surf_usm_v(1)%ncsws,                &
+                                           surf_def_v(2)%ncsws,                &
+                                           surf_lsm_v(2)%ncsws,                &
+                                           surf_usm_v(2)%ncsws,                &
+                                           surf_def_v(3)%ncsws,                &
+                                           surf_lsm_v(3)%ncsws,                &
+                                           surf_usm_v(3)%ncsws, n_out ) 
+
+            CASE ( 'qrsws' )
+               CALL surface_output_sum_up( surf_def_h(0)%qrsws,                &
+                                           surf_def_h(1)%qrsws,                &
+                                           surf_lsm_h%qrsws,                   &
+                                           surf_usm_h%qrsws,                   &
+                                           surf_def_v(0)%qrsws,                &
+                                           surf_lsm_v(0)%qrsws,                &
+                                           surf_usm_v(0)%qrsws,                &
+                                           surf_def_v(1)%qrsws,                &
+                                           surf_lsm_v(1)%qrsws,                &
+                                           surf_usm_v(1)%qrsws,                &
+                                           surf_def_v(2)%qrsws,                &
+                                           surf_lsm_v(2)%qrsws,                &
+                                           surf_usm_v(2)%qrsws,                &
+                                           surf_def_v(3)%qrsws,                &
+                                           surf_lsm_v(3)%qrsws,                &
+                                           surf_usm_v(3)%qrsws, n_out ) 
+
+            CASE ( 'nrsws' )
+               CALL surface_output_sum_up( surf_def_h(0)%nrsws,                &
+                                           surf_def_h(1)%nrsws,                &
+                                           surf_lsm_h%nrsws,                   &
+                                           surf_usm_h%nrsws,                   &
+                                           surf_def_v(0)%nrsws,                &
+                                           surf_lsm_v(0)%nrsws,                &
+                                           surf_usm_v(0)%nrsws,                &
+                                           surf_def_v(1)%nrsws,                &
+                                           surf_lsm_v(1)%nrsws,                &
+                                           surf_usm_v(1)%nrsws,                &
+                                           surf_def_v(2)%nrsws,                &
+                                           surf_lsm_v(2)%nrsws,                &
+                                           surf_usm_v(2)%nrsws,                &
+                                           surf_def_v(3)%nrsws,                &
+                                           surf_lsm_v(3)%nrsws,                &
+                                           surf_usm_v(3)%nrsws, n_out ) 
+
+            CASE ( 'sasws' )
+               CALL surface_output_sum_up( surf_def_h(0)%sasws,                &
+                                           surf_def_h(1)%sasws,                &
+                                           surf_lsm_h%sasws,                   &
+                                           surf_usm_h%sasws,                   &
+                                           surf_def_v(0)%sasws,                &
+                                           surf_lsm_v(0)%sasws,                &
+                                           surf_usm_v(0)%sasws,                &
+                                           surf_def_v(1)%sasws,                &
+                                           surf_lsm_v(1)%sasws,                &
+                                           surf_usm_v(1)%sasws,                &
+                                           surf_def_v(2)%sasws,                &
+                                           surf_lsm_v(2)%sasws,                &
+                                           surf_usm_v(2)%sasws,                &
+                                           surf_def_v(3)%sasws,                &
+                                           surf_lsm_v(3)%sasws,                &
+                                           surf_usm_v(3)%sasws, n_out ) 
+
+            CASE ( 'q_surface' )
+               CALL surface_output_sum_up( surf_def_h(0)%q_surface,            &
+                                           surf_def_h(1)%q_surface,            &
+                                           surf_lsm_h%q_surface,               &
+                                           surf_usm_h%q_surface,               &
+                                           surf_def_v(0)%q_surface,            &
+                                           surf_lsm_v(0)%q_surface,            &
+                                           surf_usm_v(0)%q_surface,            &
+                                           surf_def_v(1)%q_surface,            &
+                                           surf_lsm_v(1)%q_surface,            &
+                                           surf_usm_v(1)%q_surface,            &
+                                           surf_def_v(2)%q_surface,            &
+                                           surf_lsm_v(2)%q_surface,            &
+                                           surf_usm_v(2)%q_surface,            &
+                                           surf_def_v(3)%q_surface,            &
+                                           surf_lsm_v(3)%q_surface,            &
+                                           surf_usm_v(3)%q_surface, n_out ) 
+
+            CASE ( 'theta_surface' )
+               CALL surface_output_sum_up( surf_def_h(0)%pt_surface,           &
+                                           surf_def_h(1)%pt_surface,           &
+                                           surf_lsm_h%pt_surface,              &
+                                           surf_usm_h%pt_surface,              &
+                                           surf_def_v(0)%pt_surface,           &
+                                           surf_lsm_v(0)%pt_surface,           &
+                                           surf_usm_v(0)%pt_surface,           &
+                                           surf_def_v(1)%pt_surface,           &
+                                           surf_lsm_v(1)%pt_surface,           &
+                                           surf_usm_v(1)%pt_surface,           &
+                                           surf_def_v(2)%pt_surface,           &
+                                           surf_lsm_v(2)%pt_surface,           &
+                                           surf_usm_v(2)%pt_surface,           &
+                                           surf_def_v(3)%pt_surface,           &
+                                           surf_lsm_v(3)%pt_surface,           &
+                                           surf_usm_v(3)%pt_surface, n_out ) 
+
+            CASE ( 'thetav_surface' )
+               CALL surface_output_sum_up( surf_def_h(0)%vpt_surface,          &
+                                           surf_def_h(1)%vpt_surface,          &
+                                           surf_lsm_h%vpt_surface,             &
+                                           surf_usm_h%vpt_surface,             &
+                                           surf_def_v(0)%vpt_surface,          &
+                                           surf_lsm_v(0)%vpt_surface,          &
+                                           surf_usm_v(0)%vpt_surface,          &
+                                           surf_def_v(1)%vpt_surface,          &
+                                           surf_lsm_v(1)%vpt_surface,          &
+                                           surf_usm_v(1)%vpt_surface,          &
+                                           surf_def_v(2)%vpt_surface,          &
+                                           surf_lsm_v(2)%vpt_surface,          &
+                                           surf_usm_v(2)%vpt_surface,          &
+                                           surf_def_v(3)%vpt_surface,          &
+                                           surf_lsm_v(3)%vpt_surface,          &
+                                           surf_usm_v(3)%vpt_surface, n_out ) 
+
+            CASE ( 'rad_net' )
+               CALL surface_output_sum_up( surf_def_h(0)%rad_net,              &
+                                           surf_def_h(1)%rad_net,              &
+                                           surf_lsm_h%rad_net,                 &
+                                           surf_usm_h%rad_net,                 &
+                                           surf_def_v(0)%rad_net,              &
+                                           surf_lsm_v(0)%rad_net,              &
+                                           surf_usm_v(0)%rad_net,              &
+                                           surf_def_v(1)%rad_net,              &
+                                           surf_lsm_v(1)%rad_net,              &
+                                           surf_usm_v(1)%rad_net,              &
+                                           surf_def_v(2)%rad_net,              &
+                                           surf_lsm_v(2)%rad_net,              &
+                                           surf_usm_v(2)%rad_net,              &
+                                           surf_def_v(3)%rad_net,              &
+                                           surf_lsm_v(3)%rad_net,              &
+                                           surf_usm_v(3)%rad_net, n_out ) 
+
+            CASE ( 'rad_lw_in' )
+               CALL surface_output_sum_up( surf_def_h(0)%rad_lw_in,            &
+                                           surf_def_h(1)%rad_lw_in,            &
+                                           surf_lsm_h%rad_lw_in,               &
+                                           surf_usm_h%rad_lw_in,               &
+                                           surf_def_v(0)%rad_lw_in,            &
+                                           surf_lsm_v(0)%rad_lw_in,            &
+                                           surf_usm_v(0)%rad_lw_in,            &
+                                           surf_def_v(1)%rad_lw_in,            &
+                                           surf_lsm_v(1)%rad_lw_in,            &
+                                           surf_usm_v(1)%rad_lw_in,            &
+                                           surf_def_v(2)%rad_lw_in,            &
+                                           surf_lsm_v(2)%rad_lw_in,            &
+                                           surf_usm_v(2)%rad_lw_in,            &
+                                           surf_def_v(3)%rad_lw_in,            &
+                                           surf_lsm_v(3)%rad_lw_in,            &
+                                           surf_usm_v(3)%rad_lw_in, n_out ) 
+
+            CASE ( 'rad_lw_out' )
+               CALL surface_output_sum_up( surf_def_h(0)%rad_lw_out,           &
+                                           surf_def_h(1)%rad_lw_out,           &
+                                           surf_lsm_h%rad_lw_out,              &
+                                           surf_usm_h%rad_lw_out,              &
+                                           surf_def_v(0)%rad_lw_out,           &
+                                           surf_lsm_v(0)%rad_lw_out,           &
+                                           surf_usm_v(0)%rad_lw_out,           &
+                                           surf_def_v(1)%rad_lw_out,           &
+                                           surf_lsm_v(1)%rad_lw_out,           &
+                                           surf_usm_v(1)%rad_lw_out,           &
+                                           surf_def_v(2)%rad_lw_out,           &
+                                           surf_lsm_v(2)%rad_lw_out,           &
+                                           surf_usm_v(2)%rad_lw_out,           &
+                                           surf_def_v(3)%rad_lw_out,           &
+                                           surf_lsm_v(3)%rad_lw_out,           &
+                                           surf_usm_v(3)%rad_lw_out, n_out ) 
+
+            CASE ( 'rad_sw_in' )
+               CALL surface_output_sum_up( surf_def_h(0)%rad_sw_in,            &
+                                           surf_def_h(1)%rad_sw_in,            &
+                                           surf_lsm_h%rad_sw_in,               &
+                                           surf_usm_h%rad_sw_in,               &
+                                           surf_def_v(0)%rad_sw_in,            &
+                                           surf_lsm_v(0)%rad_sw_in,            &
+                                           surf_usm_v(0)%rad_sw_in,            &
+                                           surf_def_v(1)%rad_sw_in,            &
+                                           surf_lsm_v(1)%rad_sw_in,            &
+                                           surf_usm_v(1)%rad_sw_in,            &
+                                           surf_def_v(2)%rad_sw_in,            &
+                                           surf_lsm_v(2)%rad_sw_in,            &
+                                           surf_usm_v(2)%rad_sw_in,            &
+                                           surf_def_v(3)%rad_sw_in,            &
+                                           surf_lsm_v(3)%rad_sw_in,            &
+                                           surf_usm_v(3)%rad_sw_in, n_out ) 
+
+            CASE ( 'rad_sw_out' )
+               CALL surface_output_sum_up( surf_def_h(0)%rad_sw_out,           &
+                                           surf_def_h(1)%rad_sw_out,           &
+                                           surf_lsm_h%rad_sw_out,              &
+                                           surf_usm_h%rad_sw_out,              &
+                                           surf_def_v(0)%rad_sw_out,           &
+                                           surf_lsm_v(0)%rad_sw_out,           &
+                                           surf_usm_v(0)%rad_sw_out,           &
+                                           surf_def_v(1)%rad_sw_out,           &
+                                           surf_lsm_v(1)%rad_sw_out,           &
+                                           surf_usm_v(1)%rad_sw_out,           &
+                                           surf_def_v(2)%rad_sw_out,           &
+                                           surf_lsm_v(2)%rad_sw_out,           &
+                                           surf_usm_v(2)%rad_sw_out,           &
+                                           surf_def_v(3)%rad_sw_out,           &
+                                           surf_lsm_v(3)%rad_sw_out,           &
+                                           surf_usm_v(3)%rad_sw_out, n_out ) 
+
+            CASE ( 'ghf' )
+               CALL surface_output_sum_up( surf_def_h(0)%ghf,                  &
+                                           surf_def_h(1)%ghf,                  &
+                                           surf_lsm_h%ghf,                     &
+                                           surf_usm_h%ghf,                     &
+                                           surf_def_v(0)%ghf,                  &
+                                           surf_lsm_v(0)%ghf,                  &
+                                           surf_usm_v(0)%ghf,                  &
+                                           surf_def_v(1)%ghf,                  &
+                                           surf_lsm_v(1)%ghf,                  &
+                                           surf_usm_v(1)%ghf,                  &
+                                           surf_def_v(2)%ghf,                  &
+                                           surf_lsm_v(2)%ghf,                  &
+                                           surf_usm_v(2)%ghf,                  &
+                                           surf_def_v(3)%ghf,                  &
+                                           surf_lsm_v(3)%ghf,                  &
+                                           surf_usm_v(3)%ghf, n_out )          
+                                                                               
+            CASE ( 'r_a' )                                                     
+               CALL surface_output_sum_up( surf_def_h(0)%r_a,                  &
+                                           surf_def_h(1)%r_a,                  &
+                                           surf_lsm_h%r_a,                     &
+                                           surf_usm_h%r_a,                     &
+                                           surf_def_v(0)%r_a,                  &
+                                           surf_lsm_v(0)%r_a,                  &
+                                           surf_usm_v(0)%r_a,                  &
+                                           surf_def_v(1)%r_a,                  &
+                                           surf_lsm_v(1)%r_a,                  &
+                                           surf_usm_v(1)%r_a,                  &
+                                           surf_def_v(2)%r_a,                  &
+                                           surf_lsm_v(2)%r_a,                  &
+                                           surf_usm_v(2)%r_a,                  &
+                                           surf_def_v(3)%r_a,                  &
+                                           surf_lsm_v(3)%r_a,                  &
+                                           surf_usm_v(3)%r_a, n_out )          
+                                                                               
+            CASE ( 'r_soil' )                                                  
+               CALL surface_output_sum_up( surf_def_h(0)%r_soil,               &
+                                           surf_def_h(1)%r_soil,               &
+                                           surf_lsm_h%r_soil,                  &
+                                           surf_usm_h%r_soil,                  &
+                                           surf_def_v(0)%r_soil,               &
+                                           surf_lsm_v(0)%r_soil,               &
+                                           surf_usm_v(0)%r_soil,               &
+                                           surf_def_v(1)%r_soil,               &
+                                           surf_lsm_v(1)%r_soil,               &
+                                           surf_usm_v(1)%r_soil,               &
+                                           surf_def_v(2)%r_soil,               &
+                                           surf_lsm_v(2)%r_soil,               &
+                                           surf_usm_v(2)%r_soil,               &
+                                           surf_def_v(3)%r_soil,               &
+                                           surf_lsm_v(3)%r_soil,               &
+                                           surf_usm_v(3)%r_soil, n_out ) 
+
+            CASE ( 'r_canopy' )
+               CALL surface_output_sum_up( surf_def_h(0)%r_canopy,             &
+                                           surf_def_h(1)%r_canopy,             &
+                                           surf_lsm_h%r_canopy,                &
+                                           surf_usm_h%r_canopy,                &
+                                           surf_def_v(0)%r_canopy,             &
+                                           surf_lsm_v(0)%r_canopy,             &
+                                           surf_usm_v(0)%r_canopy,             &
+                                           surf_def_v(1)%r_canopy,             &
+                                           surf_lsm_v(1)%r_canopy,             &
+                                           surf_usm_v(1)%r_canopy,             &
+                                           surf_def_v(2)%r_canopy,             &
+                                           surf_lsm_v(2)%r_canopy,             &
+                                           surf_usm_v(2)%r_canopy,             &
+                                           surf_def_v(3)%r_canopy,             &
+                                           surf_lsm_v(3)%r_canopy,             &
+                                           surf_usm_v(3)%r_canopy, n_out ) 
+
+            CASE ( 'r_s' )
+               CALL surface_output_sum_up( surf_def_h(0)%r_s,                  &
+                                           surf_def_h(1)%r_s,                  &
+                                           surf_lsm_h%r_s,                     &
+                                           surf_usm_h%r_s,                     &
+                                           surf_def_v(0)%r_s,                  &
+                                           surf_lsm_v(0)%r_s,                  &
+                                           surf_usm_v(0)%r_s,                  &
+                                           surf_def_v(1)%r_s,                  &
+                                           surf_lsm_v(1)%r_s,                  &
+                                           surf_usm_v(1)%r_s,                  &
+                                           surf_def_v(2)%r_s,                  &
+                                           surf_lsm_v(2)%r_s,                  &
+                                           surf_usm_v(2)%r_s,                  &
+                                           surf_def_v(3)%r_s,                  &
+                                           surf_lsm_v(3)%r_s,                  &
+                                           surf_usm_v(3)%r_s, n_out ) 
+
+         END SELECT
+      ENDDO
+      
+   
+   END SUBROUTINE surface_output_averaging
+   
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Sum-up the surface data for average output variables.  
+!------------------------------------------------------------------------------!   
+   SUBROUTINE surface_output_sum_up(  var_def_h0, var_def_h1,                  &
+                                      var_lsm_h, var_usm_h,                    &
+                                      var_def_v0, var_def_v1, var_def_v2,      &
+                                      var_def_v3,                              &
+                                      var_lsm_v0, var_lsm_v1, var_lsm_v2,      &
+                                      var_lsm_v3,                              &
+                                      var_usm_v0, var_usm_v1, var_usm_v2,      &
+                                      var_usm_v3, n_out )
+   
+      IMPLICIT NONE
+
+      INTEGER(iwp) ::  m          !< running index for surface elements
+      INTEGER(iwp) ::  n_out      !< index for output variable
+      INTEGER(iwp) ::  n_surf     !< running index for surface elements
+      
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_h0 !< output variable at upward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_h1 !< output variable at downward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_h  !< output variable at upward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_h  !< output variable at upward-facing urban-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_v0 !< output variable at northward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_v1 !< output variable at southward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_v2 !< output variable at eastward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_v3 !< output variable at westward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_v0 !< output variable at northward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_v1 !< output variable at southward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_v2 !< output variable at eastward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_v3 !< output variable at westward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_v0 !< output variable at northward-facing urban-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_v1 !< output variable at southward-facing urban-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_v2 !< output variable at eastward-facing urban-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_v3 !< output variable at westward-facing urban-type surfaces
+     
+!     
+!--   Set counter variable to zero before the variable is written to 
+!--   the output array. 
+      n_surf = 0
+      
+!
+!--   Write the horizontal surfaces.
+!--   Before each the variable is written to the output data structure, first 
+!--   check if the variable for the respective surface type is defined. 
+!--   If a variable is not defined, skip the block and increment the counter
+!--   variable by the number of surface elements of this type. Usually this 
+!--   is zere, however, there might be the situation that e.g. urban surfaces
+!--   are defined but the respective variable is not allocated for this surface
+!--   type. To write the data on the exact position, increment the counter. 
+      IF ( ALLOCATED( var_def_h0 ) )  THEN
+         DO  m = 1, surf_def_h(0)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_def_h0(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_h(0)%ns  
+      ENDIF
+      IF ( ALLOCATED( var_def_h1 ) )  THEN
+         DO  m = 1, surf_def_h(1)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_def_h1(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_h(1)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_h ) )  THEN
+         DO  m = 1, surf_lsm_h%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_lsm_h(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_h%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_h ) )  THEN
+         DO  m = 1, surf_usm_h%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_usm_h(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_h%ns
+      ENDIF
+!
+!--   Write northward-facing
+      IF ( ALLOCATED( var_def_v0 ) )  THEN
+         DO  m = 1, surf_def_v(0)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_def_v0(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_v(0)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_v0 ) )  THEN
+         DO  m = 1, surf_lsm_v(0)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_lsm_v0(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_v(0)%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_v0 ) )  THEN
+         DO  m = 1, surf_usm_v(0)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_usm_v0(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_v(0)%ns
+      ENDIF
+!
+!--   Write southward-facing
+      IF ( ALLOCATED( var_def_v1 ) )  THEN
+         DO  m = 1, surf_def_v(1)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_def_v1(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_v(1)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_v1 ) )  THEN
+         DO  m = 1, surf_lsm_v(1)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_lsm_v1(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_v(1)%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_v1 ) )  THEN
+         DO  m = 1, surf_usm_v(1)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_usm_v1(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_v(1)%ns
+      ENDIF
+!
+!--   Write eastward-facing
+      IF ( ALLOCATED( var_def_v2 ) )  THEN
+         DO  m = 1, surf_def_v(2)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_def_v2(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_v(2)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_v2 ) )  THEN
+         DO  m = 1, surf_lsm_v(2)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_lsm_v2(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_v(2)%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_v2 ) )  THEN
+         DO  m = 1, surf_usm_v(2)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_usm_v2(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_v(2)%ns
+      ENDIF
+!
+!--   Write westward-facing
+      IF ( ALLOCATED( var_def_v3 ) )  THEN
+         DO  m = 1, surf_def_v(3)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_def_v3(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_v(3)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_v3 ) )  THEN
+         DO  m = 1, surf_lsm_v(3)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_lsm_v3(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_v(3)%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_v3 ) )  THEN
+         DO  m = 1, surf_usm_v(3)%ns 
+            n_surf                        = n_surf + 1
+            surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out)      &
+                                          + var_usm_v3(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_v(3)%ns
+      ENDIF
+   
+   END SUBROUTINE surface_output_sum_up
+   
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Collect the surface data from different types and different orientation.
+!------------------------------------------------------------------------------!
+   SUBROUTINE surface_output_collect( var_def_h0, var_def_h1,                  &
+                                      var_lsm_h, var_usm_h,                    &
+                                      var_def_v0, var_def_v1, var_def_v2,      &
+                                      var_def_v3,                              &
+                                      var_lsm_v0, var_lsm_v1, var_lsm_v2,      &
+                                      var_lsm_v3,                              &
+                                      var_usm_v0, var_usm_v1, var_usm_v2,      &
+                                      var_usm_v3 )
+   
+      IMPLICIT NONE
+
+      INTEGER(iwp) ::  m      !< running index for surface elements
+      INTEGER(iwp) ::  n_surf !< running index for surface elements
+      
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_h0 !< output variable at upward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_h1 !< output variable at downward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_h  !< output variable at upward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_h  !< output variable at upward-facing urban-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_v0 !< output variable at northward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_v1 !< output variable at southward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_v2 !< output variable at eastward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_def_v3 !< output variable at westward-facing default-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_v0 !< output variable at northward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_v1 !< output variable at southward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_v2 !< output variable at eastward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_lsm_v3 !< output variable at westward-facing natural-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_v0 !< output variable at northward-facing urban-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_v1 !< output variable at southward-facing urban-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_v2 !< output variable at eastward-facing urban-type surfaces
+      REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) ::  var_usm_v3 !< output variable at westward-facing urban-type surfaces
+     
+!     
+!--   Set counter variable to zero before the variable is written to 
+!--   the output array. 
+      n_surf = 0
+      
+!
+!--   Write the horizontal surfaces.
+!--   Before each the variable is written to the output data structure, first 
+!--   check if the variable for the respective surface type is defined. 
+!--   If a variable is not defined, skip the block and increment the counter
+!--   variable by the number of surface elements of this type. Usually this 
+!--   is zere, however, there might be the situation that e.g. urban surfaces
+!--   are defined but the respective variable is not allocated for this surface
+!--   type. To write the data on the exact position, increment the counter. 
+      IF ( ALLOCATED( var_def_h0 ) )  THEN
+         DO  m = 1, surf_def_h(0)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_def_h0(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_h(0)%ns  
+      ENDIF
+      IF ( ALLOCATED( var_def_h1 ) )  THEN
+         DO  m = 1, surf_def_h(1)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_def_h1(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_h(1)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_h ) )  THEN
+         DO  m = 1, surf_lsm_h%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_lsm_h(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_h%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_h ) )  THEN
+         DO  m = 1, surf_usm_h%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_usm_h(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_h%ns
+      ENDIF
+!
+!--   Write northward-facing
+      IF ( ALLOCATED( var_def_v0 ) )  THEN
+         DO  m = 1, surf_def_v(0)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_def_v0(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_v(0)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_v0 ) )  THEN
+         DO  m = 1, surf_lsm_v(0)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_lsm_v0(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_v(0)%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_v0 ) )  THEN
+         DO  m = 1, surf_usm_v(0)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_usm_v0(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_v(0)%ns
+      ENDIF
+!
+!--   Write southward-facing
+      IF ( ALLOCATED( var_def_v1 ) )  THEN
+         DO  m = 1, surf_def_v(1)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_def_v1(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_v(1)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_v1 ) )  THEN
+         DO  m = 1, surf_lsm_v(1)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_lsm_v1(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_v(1)%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_v1 ) )  THEN
+         DO  m = 1, surf_usm_v(1)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_usm_v1(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_v(1)%ns
+      ENDIF
+!
+!--   Write eastward-facing
+      IF ( ALLOCATED( var_def_v2 ) )  THEN
+         DO  m = 1, surf_def_v(2)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_def_v2(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_v(2)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_v2 ) )  THEN
+         DO  m = 1, surf_lsm_v(2)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_lsm_v2(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_v(2)%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_v2 ) )  THEN
+         DO  m = 1, surf_usm_v(2)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_usm_v2(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_v(2)%ns
+      ENDIF
+!
+!--   Write westward-facing
+      IF ( ALLOCATED( var_def_v3 ) )  THEN
+         DO  m = 1, surf_def_v(3)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_def_v3(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_def_v(3)%ns
+      ENDIF
+      IF ( ALLOCATED( var_lsm_v3 ) )  THEN
+         DO  m = 1, surf_lsm_v(3)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_lsm_v3(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_lsm_v(3)%ns
+      ENDIF
+      IF ( ALLOCATED( var_usm_v3 ) )  THEN
+         DO  m = 1, surf_usm_v(3)%ns 
+            n_surf                   = n_surf + 1
+            surfaces%var_out(n_surf) = var_usm_v3(m)
+         ENDDO
+      ELSE
+         n_surf = n_surf + surf_usm_v(3)%ns
+      ENDIF
+   
+   END SUBROUTINE surface_output_collect
+   
+   
+   
+!    SUBROUTINE data_output_surf_driver (time)
+!       IMPLICIT NONE
+! 
+!       REAL(kind=wp),INTENT(IN)    :: time
+!       CHARACTER(LEN=64)           :: filename
+!       INTEGER,SAVE                :: file_number=1
+! 
+!       if(output_mode == "legacy_VTK")   then
+!           CALL prepare_vtk_data
+! 
+!           write(filename,'(a,i5.5,a)') 'surface_output_',file_number,'.vtk'
+!           if(myid == surf_io_pe)   write(9,*) 'Legacy VTK surface output to file: ',trim(filename)
+!           file_number = file_number+1
+! 
+!           CALL write_ascii_vtk (filename, time)
+!       END IF
+! 
+!       RETURN
+!    END SUBROUTINE data_output_surf_driver
+   
+   
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Parin for output of surface parameters
+!------------------------------------------------------------------------------!
+    SUBROUTINE surface_output_parin
+
+       USE control_parameters,                                                 &
+           ONLY:  message_string
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
+       
+
+       NAMELIST /surface_output_parameters/                                    &
+                                  averaging_interval_surf, data_output_surf,   &
+                                  dt_dosurf, dt_dosurf_av,                     &
+                                  skip_time_dosurf, skip_time_dosurf_av
+                                  
+       line = ' '
+  
+!
+!--    Try to find the namelist
+       REWIND ( 11 )
+       line = ' '
+       DO WHILE ( INDEX( line, '&surface_output_parameters' ) == 0 )
+          READ ( 11, '(A)', END=14 )  line
+       ENDDO
+       BACKSPACE ( 11 )
+
+!
+!--    Read namelist
+       READ ( 11, surface_output_parameters, ERR = 10 )
+!
+!--    Set flag that indicates that surface data output is switched on
+       surface_data_output = .TRUE.
+       GOTO 14
+
+ 10    BACKSPACE( 11 )
+       READ( 11 , '(A)') line
+       CALL parin_fail_message( 'surface_output_parameters', line )      
+       
+ 14    CONTINUE
+       
+
+    END SUBROUTINE surface_output_parin
+    
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check the input parameters for consistency. Further pre-process the given 
+!> output variables, i.e. separate them into average and non-average output
+!> variables and map them onto internal output array. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE surface_output_check_parameters
+
+       USE control_parameters,                                                 &
+           ONLY:  averaging_interval, dt_data_output, dt_data_output_av,       &
+                  message_string, skip_time_data_output,                       &
+                  skip_time_data_output_av
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=100) ::  trimvar !< dummy for single output variable
+
+       INTEGER(iwp) ::  av    !< id indicating average or non-average data output
+       INTEGER(iwp) ::  ilen  !< string length
+       INTEGER(iwp) ::  n_out !< running index for number of output variables
+
+!
+!--    Check the average interval
+       IF ( averaging_interval_surf == 9999999.9_wp )  THEN
+          averaging_interval_surf = averaging_interval
+       ENDIF
+!
+!--    Set the default data-output interval dt_data_output if necessary
+       IF ( dt_dosurf    == 9999999.9_wp )  dt_dosurf    = dt_data_output
+       IF ( dt_dosurf_av == 9999999.9_wp )  dt_dosurf_av = dt_data_output
+
+       IF ( averaging_interval_surf > dt_dosurf_av )  THEN
+          WRITE( message_string, * )  'averaging_interval_surf = ',            &
+                averaging_interval_surf, ' must be <= dt_dosurf_av = ',        &
+                dt_dosurf_av
+          CALL message( 'surface_output_check_parameters',                     &
+                        'PA0087', 1, 2, 0, 6, 0 )
+       ENDIF
+
+!
+!--   Count number of output variables and separate output strings for 
+!--   average and non-average output variables.
+      n_out = 0
+      DO WHILE ( data_output_surf(n_out+1)(1:1) /= ' ' )
+
+         n_out = n_out + 1
+         ilen = LEN_TRIM( data_output_surf(n_out) )
+         trimvar = TRIM( data_output_surf(n_out) )
+
+!
+!--      Check for data averaging
+         av = 0
+         IF ( ilen > 3 )  THEN
+            IF ( data_output_surf(n_out)(ilen-2:ilen) == '_av' )  THEN
+               trimvar = data_output_surf(n_out)(1:ilen-3)
+               av      = 1
+            ENDIF
+         ENDIF
+
+         dosurf_no(av) = dosurf_no(av) + 1
+         dosurf(av,dosurf_no(av)) = TRIM( trimvar )
+
+!
+!--      Check if all output variables are known
+         SELECT CASE ( TRIM( trimvar ) )
+
+            CASE ( 'css', 'cssws', 'qsws_liq', 'qsws_soil', 'qsws_veg' )
+               message_string = TRIM( trimvar ) //                             &
+                             ' is not yet implemented in the surface output'
+               CALL message( 'surface_output_check_parameters',                &
+                             '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 ( 'usws', 'vsws' )
+
+            CASE ( 'qcsws', 'ncsws', 'qrsws', 'nrsws', 'sasws' )
+
+            CASE ( 'shf', 'qsws', 'ssws' )
+
+            CASE ( 'q_surface', 'theta_surface', 'thetav_surface' )
+
+            CASE ( 'rad_net' )
+
+            CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
+
+            CASE ( 'ghf' )
+
+            CASE ( 'r_a', 'r_canopy', 'r_soil', 'r_s' )
+
+            CASE DEFAULT
+               message_string = TRIM( trimvar ) //                             &
+                             ' is not part of the surface output'
+               CALL message( 'surface_output_check_parameters',                &
+                             'PA0087', 1, 2, 0, 6, 0 )
+         END SELECT
+
+       ENDDO
+
+    END SUBROUTINE surface_output_check_parameters
+
+!--Private Subroutines (Only called from tinside his module)
+
+
+
+!    SUBROUTINE write_ascii_vtk (filename,time)
+!       IMPLICIT NONE
+!       CHARACTER(LEN=*),INTENT(IN)   :: filename
+!       REAL(kind=wp),INTENT(IN)      :: time
+! 
+!       INTEGER                       :: file_id=251
+!       INTEGER                       :: n
+! 
+!       IF(myid == surf_io_pe)  then
+!          open(file_id,file=TRIM(filename))                                 ! Please use check_open in the final version
+! 
+!          write(file_id,'(a)') "# vtk DataFile Version 3.0"
+!          write(file_id,'(a,f8.2,a)') "legacy vtk File generated by PALM,  simulation time = ",time," sec"
+!          write(file_id,'(a)') "ASCII"
+! 
+!          write(file_id,'(a)') "DATASET POLYDATA"
+!          write(file_id,'(a,i5,a)') "POINTS ",npo," float"
+!          do n=1,npo
+!             write (file_id,'(8f10.1)') points(1:3,n)
+!          end do
+! 
+!          write(file_id,'(a,8i10)') "POLYGONS ",npg, 5*npg
+!          do n=1,npg
+!             write(file_id,'(8i10)') polygons(:,n)
+!          end do
+! 
+!          write(file_id,'(a,i10)') "CELL_DATA ",npg
+!          write(file_id,'(a,i10)') "SCALARS cell_scalars float 1 "
+!          write(file_id,'(a,i10)') "LOOKUP_TABLE default "
+! 
+!          do n=1,npg
+!             write(file_id,'(f10.2)') cell_scalars(n)
+!          end do
+! 
+!          close (file_id)
+! 
+!          DEALLOCATE (points,polygons,cell_scalars)
+!       ENDIF
+! 
+!       RETURN
+!    END SUBROUTINE write_ascii_vtk
+
+END MODULE surface_output_mod
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 3420)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Surface data output
+! 
+! 3418 2018-10-24 16:07:39Z kanani
 ! call to material_heat_model now with check if spinup runs (rvtils)
 ! 
@@ -432,6 +435,7 @@
                skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
                skip_time_do3d, skip_time_domask, skip_time_dopr,               &
-               skip_time_data_output_av, sloping_surface,                      &
-               stop_dt, terminate_coupled, terminate_run, timestep_scheme,     &
+               skip_time_data_output_av, sloping_surface, stop_dt,             &
+               surface_data_output, terminate_coupled, terminate_run,          &
+               timestep_scheme,                                                &
                time_coupling, time_do2d_xy, time_do2d_xz, time_do2d_yz,        &
                time_do3d, time_domask, time_dopr, time_dopr_av,                &
@@ -519,4 +523,10 @@
     USE surface_mod,                                                           &
         ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
+        
+    USE surface_output_mod,                                                    &
+        ONLY:  average_count_surf, averaging_interval_surf, dt_dosurf,         &
+               dt_dosurf_av, surface_output, surface_output_averaging,         &
+               skip_time_dosurf, skip_time_dosurf_av, time_dosurf,             &
+               time_dosurf_av
 
     USE turbulence_closure_mod,                                                &
@@ -1235,4 +1245,14 @@
        time_run_control   = time_run_control + dt_3d
 !
+!--    Increment time-counter for surface output
+       IF ( surface_data_output )  THEN
+          IF ( simulated_time >= skip_time_dosurf )  THEN
+             time_dosurf    = time_dosurf + dt_3d
+          ENDIF
+          IF ( simulated_time >= skip_time_dosurf_av )  THEN
+             time_dosurf_av = time_dosurf_av + dt_3d
+          ENDIF
+       ENDIF
+!
 !--    In case of synthetic turbulence generation and parametrized turbulence
 !--    information, update the time counter and if required, adjust the 
@@ -1324,4 +1344,16 @@
              average_count_3d = average_count_3d + 1
              time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
+          ENDIF
+       ENDIF
+!
+!--    Average surface data
+       IF ( surface_data_output )  THEN
+          IF ( averaging_interval_surf /= 0.0_wp  .AND.                        & 
+               ( dt_dosurf_av - time_dosurf_av ) <= averaging_interval_surf    &
+          .AND.  simulated_time >= skip_time_dosurf_av )  THEN
+             IF ( time_dosurf_av >= dt_averaging_input )  THEN        
+                CALL surface_output_averaging
+                average_count_surf = average_count_surf + 1
+             ENDIF
           ENDIF
        ENDIF
@@ -1424,4 +1456,16 @@
           time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
        ENDIF
+!
+!--    Output of surface data, instantaneous and averaged data
+       IF ( surface_data_output )  THEN
+          IF ( time_dosurf >= dt_dosurf )  THEN
+             CALL surface_output( 0 )
+             time_dosurf = MOD( time_dosurf, MAX( dt_dosurf, dt_3d ) )
+          ENDIF
+          IF ( time_dosurf_av >= dt_dosurf_av )  THEN
+             CALL surface_output( 1 )
+             time_dosurf_av = MOD( time_dosurf_av, MAX( dt_dosurf_av, dt_3d ) )
+          ENDIF
+       ENDIF
 
 !
Index: /palm/trunk/SOURCE/user_spectra.f90
===================================================================
--- /palm/trunk/SOURCE/user_spectra.f90	(revision 3420)
+++ /palm/trunk/SOURCE/user_spectra.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -105,5 +108,5 @@
        SELECT CASE ( TRIM( data_output_sp(m) ) )
           
-          CASE ( 'u', 'v', 'w', 'pt', 'q', 's' )
+          CASE ( 'u', 'v', 'w', 'theta', 'q', 's' )
 !--          Not allowed here since these are the standard quantities used in 
 !--          preprocess_spectra.
@@ -124,5 +127,5 @@
        SELECT CASE ( TRIM( data_output_sp(m) ) )
 
-          CASE ( 'u', 'v', 'w', 'pt', 'q', 's' )
+          CASE ( 'u', 'v', 'w', 'theta', 'q', 's' )
 !--          Not allowed here since these are the standard quantities used in 
 !--          data_output_spectra.
Index: /palm/trunk/SOURCE/virtual_flight_mod.f90
===================================================================
--- /palm/trunk/SOURCE/virtual_flight_mod.f90	(revision 3420)
+++ /palm/trunk/SOURCE/virtual_flight_mod.f90	(revision 3421)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed output variables
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -496,5 +499,5 @@
                 k             = k + 1
              ELSEIF ( i == id_pt )  THEN       
-                dofl_label(k) = TRIM( label_leg ) // '_pt'
+                dofl_label(k) = TRIM( label_leg ) // '_theta'
                 dofl_unit(k)  = 'K'
                 k             = k + 1
Index: /palm/trunk/TESTS/cases/example_cbl/INPUT/example_cbl_p3d
===================================================================
--- /palm/trunk/TESTS/cases/example_cbl/INPUT/example_cbl_p3d	(revision 3420)
+++ /palm/trunk/TESTS/cases/example_cbl/INPUT/example_cbl_p3d	(revision 3421)
@@ -74,13 +74,13 @@
                                  'w_xz',
                                  'w_xz_av',
-                                 'pt_xy',
-                                 'pt_xz',
+                                 'theta_xy',
+                                 'theta_xz',
 
-    data_output_pr             = 'wpt',
-                                 'w"pt"',
-                                 'w*pt*',
+    data_output_pr             = 'wtheta',
+                                 'w"theta"',
+                                 'w*theta*',
                                  'w*2',
-                                 '#pt',
-                                 'pt*2',
+                                 '#theta',
+                                 'theta*2',
 
     section_xy                 = 2,
Index: /palm/trunk/TESTS/cases/example_cbl_restart/INPUT/example_cbl_restart_p3dr
===================================================================
--- /palm/trunk/TESTS/cases/example_cbl_restart/INPUT/example_cbl_restart_p3dr	(revision 3420)
+++ /palm/trunk/TESTS/cases/example_cbl_restart/INPUT/example_cbl_restart_p3dr	(revision 3421)
@@ -74,13 +74,13 @@
                                  'w_xz',
                                  'w_xz_av',
-                                 'pt_xy',
-                                 'pt_xz',
+                                 'theta_xy',
+                                 'theta_xz',
 
-    data_output_pr             = 'wpt',
-                                 'w"pt"',
-                                 'w*pt*',
+    data_output_pr             = 'wtheta',
+                                 'w"theta"',
+                                 'w*theta*',
                                  'w*2',
-                                 '#pt',
-                                 'pt*2',
+                                 '#theta',
+                                 'theta*2',
 
     section_xy                 = 2,
Index: /palm/trunk/TESTS/cases/example_cbl_restart2/INPUT/example_cbl_restart2_p3dr
===================================================================
--- /palm/trunk/TESTS/cases/example_cbl_restart2/INPUT/example_cbl_restart2_p3dr	(revision 3420)
+++ /palm/trunk/TESTS/cases/example_cbl_restart2/INPUT/example_cbl_restart2_p3dr	(revision 3421)
@@ -74,13 +74,13 @@
                                  'w_xz',
                                  'w_xz_av',
-                                 'pt_xy',
-                                 'pt_xz',
+                                 'theta_xy',
+                                 'theta_xz',
 
-    data_output_pr             = 'wpt',
-                                 'w"pt"',
-                                 'w*pt*',
+    data_output_pr             = 'wtheta',
+                                 'w"theta"',
+                                 'w*theta*',
                                  'w*2',
-                                 '#pt',
-                                 'pt*2',
+                                 '#theta',
+                                 'theta*2',
 
     section_xy                 = 2,
Index: /palm/trunk/TESTS/cases/example_cbl_short/INPUT/example_cbl_short_p3d
===================================================================
--- /palm/trunk/TESTS/cases/example_cbl_short/INPUT/example_cbl_short_p3d	(revision 3420)
+++ /palm/trunk/TESTS/cases/example_cbl_short/INPUT/example_cbl_short_p3d	(revision 3421)
@@ -74,10 +74,10 @@
                                  'w_xz',
                                  'w_xz_av',
-                                 'pt_xy',
-                                 'pt_xz',
+                                 'theta_xy',
+                                 'theta_xz',
 
-    data_output_pr             = 'wpt', 'w"pt"', 'w*pt*',
+    data_output_pr             = 'wtheta', 'w"theta"', 'w*theta*',
                                  'w*2',
-                                 '#pt', 'pt*2',
+                                 '#theta', 'theta*2',
 
     section_xy                 = 2,
Index: /palm/trunk/TESTS/cases/example_cbl_user_code/INPUT/example_cbl_user_code_p3d
===================================================================
--- /palm/trunk/TESTS/cases/example_cbl_user_code/INPUT/example_cbl_user_code_p3d	(revision 3420)
+++ /palm/trunk/TESTS/cases/example_cbl_user_code/INPUT/example_cbl_user_code_p3d	(revision 3421)
@@ -74,13 +74,13 @@
                                  'w_xz',
                                  'w_xz_av',
-                                 'pt_xy',
-                                 'pt_xz',
+                                 'theta_xy',
+                                 'theta_xz',
 
-    data_output_pr             = 'wpt',
-                                 'w"pt"',
-                                 'w*pt*',
+    data_output_pr             = 'wtheta',
+                                 'w"theta"',
+                                 'w*theta*',
                                  'w*2',
-                                 '#pt',
-                                 'pt*2',
+                                 '#theta',
+                                 'theta*2',
 
     section_xy                 = 2,
Index: /palm/trunk/TESTS/cases/urban_environment/INPUT/urban_environment_p3d
===================================================================
--- /palm/trunk/TESTS/cases/urban_environment/INPUT/urban_environment_p3d	(revision 3420)
+++ /palm/trunk/TESTS/cases/urban_environment/INPUT/urban_environment_p3d	(revision 3421)
@@ -91,5 +91,5 @@
                             '#v',  'v*2',  'wv',  'w*v*',  'w"v"',
                             'w',   'w*2',
-                            '#pt', 'pt*2', 'wpt', 'w*pt*', 'w"pt"',
+                            '#theta', 'theta*2', 'wtheta', 'w*theta*', 'w"theta"',
                             '#q',  'q*2',
                             'e', 'e*', '#km', '#l',
@@ -111,5 +111,5 @@
                   'v',  'v_av',
                   'w',  'w_av',
-                  'pt', 'pt_av',
+                  'theta', 'theta_av',
                   'q',  'q_av',
                   'e',  'e_av',
@@ -119,5 +119,5 @@
                   'shf*_xy', 'shf*_xy_av',
                   'ghf*_xy', 'ghf*_xy_av',
-                  'u*_xy', 'z0*_xy',
+                  'us*_xy', 'z0*_xy',
                   't*_xy', 'ol*_xy',
                   'm_liq*_xy',
Index: /palm/trunk/TESTS/cases/urban_environment_restart/INPUT/urban_environment_restart_p3dr
===================================================================
--- /palm/trunk/TESTS/cases/urban_environment_restart/INPUT/urban_environment_restart_p3dr	(revision 3420)
+++ /palm/trunk/TESTS/cases/urban_environment_restart/INPUT/urban_environment_restart_p3dr	(revision 3421)
@@ -91,5 +91,5 @@
                             '#v',  'v*2',  'wv',  'w*v*',  'w"v"',
                             'w',   'w*2',
-                            '#pt', 'pt*2', 'wpt', 'w*pt*', 'w"pt"',
+                            '#theta', 'theta*2', 'wtheta', 'w*theta*', 'w"theta"',
                             '#q',  'q*2',
                             'e', 'e*', '#km', '#l',
@@ -111,5 +111,5 @@
                   'v',  'v_av',
                   'w',  'w_av',
-                  'pt', 'pt_av',
+                  'theta', 'theta_av',
                   'q',  'q_av',
                   'e',  'e_av',
@@ -119,5 +119,5 @@
                   'shf*_xy', 'shf*_xy_av',
                   'ghf*_xy', 'ghf*_xy_av',
-                  'u*_xy', 'z0*_xy',
+                  'us*_xy', 'z0*_xy',
                   't*_xy', 'ol*_xy',
                   'm_liq*_xy',
