Index: palm/trunk/SOURCE/surface_data_output_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_data_output_mod.f90	(revision 3735)
+++ palm/trunk/SOURCE/surface_data_output_mod.f90	(revision 3736)
@@ -25,41 +25,45 @@
 ! -----------------
 ! $Id$
+! Add azimuth and zenith to output file; set long-name attributes; 
+! clean-up coding layout
+! 
+! 3735 2019-02-12 09:52:40Z dom_dwd_user
 ! - Split initialization into initialization of arrays and further initialization
 !   in order to enable reading of restart data.
 ! - Consider restarts in surface data averaging.
 ! - Correct error message numbers
-! 
+!
 ! 3731 2019-02-11 13:06:27Z suehring
 ! Bugfix: add cpp options
-! 
+!
 ! 3727 2019-02-08 14:52:10Z gronemeier
 ! Enable NetCDF output for surface data (suehring, gronemeier)
-! 
+!
 ! 3691 2019-01-23 09:57:04Z suehring
-! Add output of surface-parallel flow speed 
-! 
+! Add output of surface-parallel flow speed
+!
 ! 3648 2019-01-02 16:35:46Z suehring
 ! Rename module and subroutines
-! 
+!
 ! 3646 2018-12-28 17:58:49Z kanani
 ! Bugfix: use time_since_reference_point instead of simulated_time (relevant
 ! when using wall/soil spinup)
-! 
+!
 ! 3614 2018-12-10 07:05:46Z raasch
 ! unused variables removed
-! 
+!
 ! 3572 2018-11-28 11:40:28Z suehring
 ! Added short- and longwave radiation flux arrays (e.g. diffuse, direct,
 ! reflected, resedual) for all surfaces (M. Salim)
-! 
+!
 ! 3494 2018-11-06 14:51:27Z suehring
-! Bugfix in gathering surface data from different types and orientation. 
+! Bugfix in gathering surface data from different types and orientation.
 ! Output of total number of surfaces and vertices added.
 ! String length is output, for more convinient post-processing.
 ! Last actions added.
-! 
+!
 ! 3483 2018-11-02 14:19:26Z raasch
 ! bugfix: missed directives for MPI added
-! 
+!
 ! 3421 2018-10-24 18:39:32Z gronemeier
 ! Add output variables
@@ -77,6 +81,6 @@
 ! ------------
 !> Generate output for surface data.
-!> 
-!> @todo Create namelist file for post-processing tool. 
+!>
+!> @todo Create namelist file for post-processing tool.
 !------------------------------------------------------------------------------!
 
@@ -87,16 +91,16 @@
    USE arrays_3d,                                                              &
        ONLY:  zu, zw
-       
+
    USE control_parameters,                                                     &
        ONLY:  coupling_char, data_output_during_spinup, end_time,              &
               message_string, run_description_header, simulated_time_at_begin, &
               spinup_time, surface_output
-       
+
    USE grid_variables,                                                         &
        ONLY: dx,dy
-       
+
    USE indices,                                                                &
        ONLY: nxl, nxr, nys, nyn, nzb, nzt
-       
+
 #if defined( __netcdf )
    USE NETCDF
@@ -109,7 +113,7 @@
        ONLY:  netcdf_create_att, netcdf_create_dim, netcdf_create_file,        &
               netcdf_create_global_atts, netcdf_create_var, netcdf_handle_error
-       
+
    USE pegrid
-       
+
    USE surface_mod,                                                            &
        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,                  &
@@ -119,21 +123,21 @@
 
    TYPE surf_out                      !< data structure which contains all surfaces elements of all types on subdomain
-   
+
       INTEGER(iwp) ::  ns             !< number of surface elements on subdomain
       INTEGER(iwp) ::  ns_total       !< total number of surface elements
       INTEGER(iwp) ::  npoints        !< number of points / vertices which define a surface element (on subdomain)
       INTEGER(iwp) ::  npoints_total  !< total number of points / vertices which define a surface element
-      
+
       INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  s        !< coordinate for NetCDF output, number of the surface element
-      
-      REAL(wp) ::  fillvalue = -9999.9_wp !< fillvalue for surface elements which are not defined
-      
-      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  azimuth  !< azimuth orientation coordinate for NetCDF output 
-      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  es_utm   !< E-UTM coordinate for NetCDF output 
-      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  ns_utm   !< E-UTM coordinate for NetCDF output 
-      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  xs       !< x-coordinate for NetCDF output 
-      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  ys       !< y-coordinate for NetCDF output 
-      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zs       !< z-coordinate for NetCDF output 
-      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zenith   !< zenith orientation coordinate for NetCDF output 
+
+      REAL(wp) ::  fillvalue = -9999.0_wp !< fillvalue for surface elements which are not defined
+
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  azimuth  !< azimuth orientation coordinate for NetCDF output
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  es_utm   !< E-UTM coordinate for NetCDF output
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  ns_utm   !< E-UTM coordinate for NetCDF output
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  xs       !< x-coordinate for NetCDF output
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  ys       !< y-coordinate for NetCDF output
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zs       !< z-coordinate for NetCDF output
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zenith   !< zenith orientation coordinate for NetCDF output
       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_out  !< output variables
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_av   !< variables used for averaging
@@ -141,19 +145,20 @@
       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 
+   CHARACTER(LEN=100), DIMENSION(0:1,300)   ::  dosurf = ' '            !< internal variable which describes the output variables
                                                                         !<  and separates averaged from non-averaged output
    CHARACTER(LEN=100), DIMENSION(0:1,300)   ::  dosurf_unit = ' '       !< internal variable which holds the unit of the given output variable
-   
+
    INTEGER(iwp) ::  average_count_surf = 0   !< number of ensemble members used for averaging
    INTEGER(iwp) ::  dosurf_no(0:1) = 0       !< number of surface output quantities
 
    INTEGER(iwp) :: nc_stat                   !< error code for netcdf routines
-   INTEGER(iwp), SAVE ::  oldmode            !< save old set-fill-mode of netcdf file (not needed, but required for routine call)
+   INTEGER(iwp) ::  oldmode                  !< save old set-fill-mode of netcdf file (not needed, but required for routine call)
 
    INTEGER(iwp), DIMENSION(0:1) ::  id_dim_s_surf         !< netcdf ID for dimension s
    INTEGER(iwp), DIMENSION(0:1) ::  id_dim_time_surf      !< netcdf ID for dimension time
    INTEGER(iwp), DIMENSION(0:1) ::  id_set_surf           !< netcdf ID for file
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_azimuth_surf   !< netcdf ID for variable azimuth
    INTEGER(iwp), DIMENSION(0:1) ::  id_var_etum_surf      !< netcdf ID for variable Es_UTM
    INTEGER(iwp), DIMENSION(0:1) ::  id_var_nutm_surf      !< netcdf ID for variable Ns_UTM
@@ -162,13 +167,14 @@
    INTEGER(iwp), DIMENSION(0:1) ::  id_var_xs_surf        !< netcdf ID for variable xs
    INTEGER(iwp), DIMENSION(0:1) ::  id_var_ys_surf        !< netcdf ID for variable ys
+   INTEGER(iwp), DIMENSION(0:1) ::  id_var_zenith_surf    !< netcdf ID for variable zenith
    INTEGER(iwp), DIMENSION(0:1) ::  id_var_zs_surf        !< netcdf ID for variable zs
    INTEGER(iwp), DIMENSION(0:1) ::  dosurf_time_count = 0 !< count of output time steps
    INTEGER(iwp), DIMENSION(0:1) ::  ntdim_surf            !< number of output time steps
-   
+
    INTEGER(iwp), DIMENSION(0:1,300) ::  id_var_dosurf     !< netcdf ID for output variables
 
    LOGICAL :: first_output(0:1) = .FALSE.                 !< true if first output was already called
-   LOGICAL :: to_netcdf = .TRUE.                          !< flag indicating parallel NetCDF output 
-   LOGICAL :: to_vtk = .TRUE.                             !< flag indicating binary surface-data output that can be further 
+   LOGICAL :: to_netcdf = .TRUE.                          !< flag indicating parallel NetCDF output
+   LOGICAL :: to_vtk = .TRUE.                             !< flag indicating binary surface-data output that can be further
                                                           !< processed to VTK format
 
@@ -182,5 +188,5 @@
 
    TYPE(surf_out) ::  surfaces      !< variable which contains all required output information
-   
+
    SAVE
 
@@ -190,41 +196,41 @@
       MODULE PROCEDURE surface_data_output
    END INTERFACE  surface_data_output
-   
+
    INTERFACE  surface_data_output_averaging
       MODULE PROCEDURE surface_data_output_averaging
    END INTERFACE  surface_data_output_averaging
-   
+
    INTERFACE  surface_data_output_check_parameters
       MODULE PROCEDURE surface_data_output_check_parameters
    END INTERFACE  surface_data_output_check_parameters
-   
+
    INTERFACE  surface_data_output_init
       MODULE PROCEDURE surface_data_output_init
    END INTERFACE  surface_data_output_init
-   
+
    INTERFACE  surface_data_output_init_arrays
       MODULE PROCEDURE surface_data_output_init_arrays
    END INTERFACE  surface_data_output_init_arrays
-   
+
    INTERFACE  surface_data_output_last_action
       MODULE PROCEDURE surface_data_output_last_action
    END INTERFACE  surface_data_output_last_action
-      
+
    INTERFACE  surface_data_output_parin
       MODULE PROCEDURE surface_data_output_parin
    END INTERFACE  surface_data_output_parin
-   
+
    INTERFACE  surface_data_output_rrd_global
       MODULE PROCEDURE surface_data_output_rrd_global
    END INTERFACE  surface_data_output_rrd_global
-   
+
    INTERFACE  surface_data_output_rrd_local
       MODULE PROCEDURE surface_data_output_rrd_local
    END INTERFACE  surface_data_output_rrd_local
-   
+
    INTERFACE  surface_data_output_wrd_global
       MODULE PROCEDURE surface_data_output_wrd_global
    END INTERFACE  surface_data_output_wrd_global
-   
+
    INTERFACE  surface_data_output_wrd_local
       MODULE PROCEDURE surface_data_output_wrd_local
@@ -250,8 +256,8 @@
 ! ------------
 !> This routine counts the number of surfaces on each core and allocates
-!> arrays. 
+!> arrays.
 !------------------------------------------------------------------------------!
    SUBROUTINE surface_data_output_init_arrays
-   
+
       IMPLICIT NONE
 
@@ -259,9 +265,9 @@
 !--   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_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   
+                  + 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
 !
 !--    Determine the total number of surfaces in the model domain
@@ -277,5 +283,5 @@
       surfaces%var_out = surfaces%fillvalue
 !
-!--   If there is an output of time average output variables, allocate the 
+!--   If there is an output of time average output variables, allocate the
 !--   required array.
       IF ( dosurf_no(1) > 0 )  THEN
@@ -283,8 +289,8 @@
          surfaces%var_av = 0.0_wp
       ENDIF
-      
+
    END SUBROUTINE surface_data_output_init_arrays
- 
- 
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -294,5 +300,5 @@
 !------------------------------------------------------------------------------!
    SUBROUTINE surface_data_output_init
-   
+
       IMPLICIT NONE
 
@@ -311,9 +317,9 @@
       INTEGER(iwp) ::  point_index_count !< local counter variable for point index
       INTEGER(iwp) ::  start_count       !< local start counter for the surface index
-      
+
       INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_points_on_pe   !< array which contains the number of points on all mpi ranks
       INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_surfaces_on_pe !< array which contains the number of surfaces on all mpi ranks
       INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) ::  point_index !< dummy array used to check where the reference points for surface polygons are located
-           
+
       REAL(wp) ::  az    !< azimuth angle, indicated the vertical orientation of a surface element
       REAL(wp) ::  off_x !< grid offset in x-direction between the stored grid index and the actual wall
@@ -323,6 +329,6 @@
 
 !
-!--   If output to VTK format is enabled, initialize point and polygon data. 
-!--   In a first step, count the number of points which are defining 
+!--   If output to VTK format is enabled, initialize point and polygon data.
+!--   In a first step, count the number of points which are defining
 !--   the surfaces and the polygons.
       IF ( to_vtk )  THEN
@@ -340,22 +346,22 @@
                k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
 !
-!--            Check if the vertices that define the surface element are already 
+!--            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    
+                  surfaces%npoints   = surfaces%npoints + 1
                   point_index(k,j,i) = surfaces%npoints - 1
                ENDIF
                IF ( point_index(k,j,i+1) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k,j,i+1) = surfaces%npoints - 1
                ENDIF
                IF ( point_index(k,j+1,i+1) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
+                  surfaces%npoints       = surfaces%npoints + 1
                   point_index(k,j+1,i+1) = surfaces%npoints - 1
                ENDIF
                IF ( point_index(k,j+1,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k,j+1,i) = surfaces%npoints - 1
-               ENDIF           
+               ENDIF
             ENDDO
          ENDDO
@@ -366,19 +372,19 @@
 
             IF ( point_index(k,j,i) < 0 )  THEN
-               surfaces%npoints   = surfaces%npoints + 1    
+               surfaces%npoints   = surfaces%npoints + 1
                point_index(k,j,i) = surfaces%npoints - 1
             ENDIF
             IF ( point_index(k,j,i+1) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
+               surfaces%npoints     = surfaces%npoints + 1
                point_index(k,j,i+1) = surfaces%npoints - 1
             ENDIF
             IF ( point_index(k,j+1,i+1) < 0 )  THEN
-               surfaces%npoints       = surfaces%npoints + 1    
+               surfaces%npoints       = surfaces%npoints + 1
                point_index(k,j+1,i+1) = surfaces%npoints - 1
             ENDIF
             IF ( point_index(k,j+1,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
+               surfaces%npoints     = surfaces%npoints + 1
                point_index(k,j+1,i) = surfaces%npoints - 1
-            ENDIF    
+            ENDIF
          ENDDO
          DO  m = 1, surf_usm_h%ns
@@ -388,19 +394,19 @@
 
             IF ( point_index(k,j,i) < 0 )  THEN
-               surfaces%npoints   = surfaces%npoints + 1    
+               surfaces%npoints   = surfaces%npoints + 1
                point_index(k,j,i) = surfaces%npoints - 1
             ENDIF
             IF ( point_index(k,j,i+1) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
+               surfaces%npoints     = surfaces%npoints + 1
                point_index(k,j,i+1) = surfaces%npoints - 1
             ENDIF
             IF ( point_index(k,j+1,i+1) < 0 )  THEN
-               surfaces%npoints       = surfaces%npoints + 1    
+               surfaces%npoints       = surfaces%npoints + 1
                point_index(k,j+1,i+1) = surfaces%npoints - 1
             ENDIF
             IF ( point_index(k,j+1,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
+               surfaces%npoints     = surfaces%npoints + 1
                point_index(k,j+1,i) = surfaces%npoints - 1
-            ENDIF     
+            ENDIF
          ENDDO
 !
@@ -410,6 +416,6 @@
 !
 !--            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  
+!--            topography. Please note, j-index for north-facing surfaces
+!--            ( l==0 ) is identical to the reference j-index outside the grid
 !--            box. Equivalent for east-facing surfaces and i-index.
                i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
@@ -419,16 +425,16 @@
 !--            Lower left /front vertex
                IF ( point_index(k,j,i) < 0 )  THEN
-                  surfaces%npoints   = surfaces%npoints + 1    
+                  surfaces%npoints   = surfaces%npoints + 1
                   point_index(k,j,i) = surfaces%npoints - 1
-               ENDIF 
+               ENDIF
 !
 !--            Upper / lower right index for north- and south-facing surfaces
                IF ( l == 0  .OR.  l == 1 )  THEN
                   IF ( point_index(k,j,i+1) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j,i+1) = surfaces%npoints - 1
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j,i+1) = surfaces%npoints - 1
                   ENDIF
@@ -437,9 +443,9 @@
                ELSEIF ( l == 2  .OR.  l == 3 )  THEN
                   IF ( point_index(k,j+1,i) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j+1,i) = surfaces%npoints - 1
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j+1,i) = surfaces%npoints - 1
                   ENDIF
@@ -448,5 +454,5 @@
 !--            Upper left / front vertex
                IF ( point_index(k+1,j,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k+1,j,i) = surfaces%npoints - 1
                ENDIF
@@ -459,16 +465,16 @@
 !--            Lower left /front vertex
                IF ( point_index(k,j,i) < 0 )  THEN
-                  surfaces%npoints   = surfaces%npoints + 1    
+                  surfaces%npoints   = surfaces%npoints + 1
                   point_index(k,j,i) = surfaces%npoints - 1
-               ENDIF 
+               ENDIF
 !
 !--            Upper / lower right index for north- and south-facing surfaces
                IF ( l == 0  .OR.  l == 1 )  THEN
                   IF ( point_index(k,j,i+1) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j,i+1) = surfaces%npoints - 1
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j,i+1) = surfaces%npoints - 1
                   ENDIF
@@ -477,9 +483,9 @@
                ELSEIF ( l == 2  .OR.  l == 3 )  THEN
                   IF ( point_index(k,j+1,i) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j+1,i) = surfaces%npoints - 1
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j+1,i) = surfaces%npoints - 1
                   ENDIF
@@ -488,5 +494,5 @@
 !--            Upper left / front vertex
                IF ( point_index(k+1,j,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k+1,j,i) = surfaces%npoints - 1
                ENDIF
@@ -500,16 +506,16 @@
 !--            Lower left /front vertex
                IF ( point_index(k,j,i) < 0 )  THEN
-                  surfaces%npoints   = surfaces%npoints + 1    
+                  surfaces%npoints   = surfaces%npoints + 1
                   point_index(k,j,i) = surfaces%npoints - 1
-               ENDIF 
+               ENDIF
 !
 !--            Upper / lower right index for north- and south-facing surfaces
                IF ( l == 0  .OR.  l == 1 )  THEN
                   IF ( point_index(k,j,i+1) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j,i+1) = surfaces%npoints - 1
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j,i+1) = surfaces%npoints - 1
                   ENDIF
@@ -518,9 +524,9 @@
                ELSEIF ( l == 2  .OR.  l == 3 )  THEN
                   IF ( point_index(k,j+1,i) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j+1,i) = surfaces%npoints - 1
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j+1,i) = surfaces%npoints - 1
                   ENDIF
@@ -529,5 +535,5 @@
 !--            Upper left / front vertex
                IF ( point_index(k+1,j,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k+1,j,i) = surfaces%npoints - 1
                ENDIF
@@ -536,14 +542,14 @@
          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 the number of points and polygons. Note, the number of
+!--      polygons is identical to the number of surfaces elements, whereas the
+!--      number of points (vertices), which define the polygons, can be larger.
          ALLOCATE( surfaces%points(3,1:surfaces%npoints) )
          ALLOCATE( surfaces%polygons(5,1:surfaces%ns)    )
 !
 !--      Note, PARAVIEW expects consecutively ordered points, in order to
-!--      unambiguously identify surfaces. 
+!--      unambiguously identify surfaces.
 !--      Hence, all PEs should know where they start counting, depending on the
-!--      number of points on the other PE's with lower MPI rank. 
+!--      number of points on the other PE's with lower MPI rank.
 #if defined( __parallel )
          CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER,                 &
@@ -554,6 +560,6 @@
 
 !
-!--      After the number of vertices is counted, repeat the loops and define the
-!--      vertices. Start with the horizontal default surfaces.
+!--      After the number of vertices is counted, repeat the loops and define
+!--      the vertices. Start with the horizontal default surfaces.
 !--      First, however, determine the offset where couting of points should be
 !--      started, which is the sum of points of all PE's with lower MPI rank.
@@ -564,13 +570,13 @@
             i                 = i + 1
          ENDDO
-          
+
          surfaces%npoints = 0
          point_index      = -1
          npg              = 0
-      
+
          DO  l = 0, 1
             DO  m = 1, surf_def_h(0)%ns
 !
-!--            Determine the indices of the respective grid cell inside the 
+!--            Determine the indices of the respective grid cell inside the
 !--            topography.
                i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
@@ -578,10 +584,10 @@
                k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
 !
-!--            Check if the vertices that define the surface element are  
+!--            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    
+                  surfaces%npoints   = surfaces%npoints + 1
                   point_index(k,j,i) = point_index_count
-                  point_index_count  = point_index_count + 1               
+                  point_index_count  = point_index_count + 1
                   surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                   surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
@@ -589,7 +595,7 @@
                ENDIF
                IF ( point_index(k,j,i+1) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k,j,i+1) = point_index_count
-                  point_index_count    = point_index_count + 1 
+                  point_index_count    = point_index_count + 1
                   surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
                   surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
@@ -597,7 +603,7 @@
                ENDIF
                IF ( point_index(k,j+1,i+1) < 0 )  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
+                  surfaces%npoints       = surfaces%npoints + 1
                   point_index(k,j+1,i+1) = point_index_count
-                  point_index_count      = point_index_count + 1 
+                  point_index_count      = point_index_count + 1
                   surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
                   surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
@@ -605,12 +611,12 @@
                ENDIF
                IF ( point_index(k,j+1,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k,j+1,i) = point_index_count
-                  point_index_count    = point_index_count + 1 
+                  point_index_count    = point_index_count + 1
                   surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
                   surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k)
                ENDIF
-               
+
                npg                        = npg + 1
                surfaces%polygons(1,npg)   = 4
@@ -626,7 +632,7 @@
             k = surf_lsm_h%k(m) + surf_lsm_h%koff
             IF ( point_index(k,j,i) < 0 )  THEN
-               surfaces%npoints   = surfaces%npoints + 1    
+               surfaces%npoints   = surfaces%npoints + 1
                point_index(k,j,i) = point_index_count
-               point_index_count  = point_index_count + 1 
+               point_index_count  = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
@@ -634,7 +640,7 @@
             ENDIF
             IF ( point_index(k,j,i+1) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
+               surfaces%npoints     = surfaces%npoints + 1
                point_index(k,j,i+1) = point_index_count
-               point_index_count    = point_index_count + 1 
+               point_index_count    = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
@@ -642,7 +648,7 @@
             ENDIF
             IF ( point_index(k,j+1,i+1) < 0 )  THEN
-               surfaces%npoints       = surfaces%npoints + 1    
+               surfaces%npoints       = surfaces%npoints + 1
                point_index(k,j+1,i+1) = point_index_count
-               point_index_count      = point_index_count + 1 
+               point_index_count      = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
@@ -650,12 +656,12 @@
             ENDIF
             IF ( point_index(k,j+1,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
+               surfaces%npoints     = surfaces%npoints + 1
                point_index(k,j+1,i) = point_index_count
-               point_index_count    = point_index_count + 1 
+               point_index_count    = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
                surfaces%points(3,surfaces%npoints) = zw(k)
             ENDIF
-            
+
             npg                        = npg + 1
             surfaces%polygons(1,npg)   = 4
@@ -663,16 +669,16 @@
             surfaces%polygons(3,npg)   = point_index(k,j,i+1)
             surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
-            surfaces%polygons(5,npg)   = point_index(k,j+1,i) 
+            surfaces%polygons(5,npg)   = point_index(k,j+1,i)
          ENDDO
-     
+
          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    
+               surfaces%npoints   = surfaces%npoints + 1
                point_index(k,j,i) = point_index_count
-               point_index_count  = point_index_count + 1 
+               point_index_count  = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
@@ -680,7 +686,7 @@
             ENDIF
             IF ( point_index(k,j,i+1) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
+               surfaces%npoints     = surfaces%npoints + 1
                point_index(k,j,i+1) = point_index_count
-               point_index_count    = point_index_count + 1 
+               point_index_count    = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
@@ -688,7 +694,7 @@
             ENDIF
             IF ( point_index(k,j+1,i+1) < 0 )  THEN
-               surfaces%npoints       = surfaces%npoints + 1    
+               surfaces%npoints       = surfaces%npoints + 1
                point_index(k,j+1,i+1) = point_index_count
-               point_index_count      = point_index_count + 1 
+               point_index_count      = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
@@ -696,12 +702,12 @@
             ENDIF
             IF ( point_index(k,j+1,i) < 0 )  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
+               surfaces%npoints     = surfaces%npoints + 1
                point_index(k,j+1,i) = point_index_count
-               point_index_count    = point_index_count + 1 
+               point_index_count    = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
                surfaces%points(3,surfaces%npoints) = zw(k)
             ENDIF
-            
+
             npg                        = npg + 1
             surfaces%polygons(1,npg)   = 4
@@ -709,32 +715,33 @@
             surfaces%polygons(3,npg)   = point_index(k,j,i+1)
             surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
-            surfaces%polygons(5,npg)   = point_index(k,j+1,i)    
+            surfaces%polygons(5,npg)   = point_index(k,j+1,i)
          ENDDO
 
          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. 
+!
+!--            Determine the indices of the respective grid cell inside the
+!--            topography.
+!--            NOTE, j-index for north-facing surfaces ( l==0 ) is
+!--            identical to the reference j-index outside the grid box.
 !--            Equivalent for east-facing surfaces and i-index.
                i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
                j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
                k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
-!        
+!
 !--            Lower left /front vertex
                IF ( point_index(k,j,i) < 0 )  THEN
-                  surfaces%npoints   = surfaces%npoints + 1    
+                  surfaces%npoints   = surfaces%npoints + 1
                   point_index(k,j,i) = point_index_count
-                  point_index_count  = point_index_count + 1 
+                  point_index_count  = point_index_count + 1
                   surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                   surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k-1)
-               ENDIF 
-!        
+               ENDIF
+!
 !--            Upper / lower right index for north- and south-facing surfaces
                IF ( l == 0  .OR.  l == 1 )  THEN
                   IF ( point_index(k,j,i+1) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j,i+1) = point_index_count
                      point_index_count    = point_index_count + 1
@@ -742,7 +749,7 @@
                      surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
                      surfaces%points(3,surfaces%npoints) = zw(k-1)
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j,i+1) = point_index_count
                      point_index_count      = point_index_count + 1
@@ -751,9 +758,9 @@
                      surfaces%points(3,surfaces%npoints) = zw(k)
                   ENDIF
-!        
+!
 !--            Upper / lower front index for east- and west-facing surfaces
                ELSEIF ( l == 2  .OR.  l == 3 )  THEN
                   IF ( point_index(k,j+1,i) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j+1,i) = point_index_count
                      point_index_count    = point_index_count + 1
@@ -761,7 +768,7 @@
                      surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
                      surfaces%points(3,surfaces%npoints) = zw(k-1)
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j+1,i) = point_index_count
                      point_index_count      = point_index_count + 1
@@ -771,8 +778,8 @@
                   ENDIF
                ENDIF
-!        
+!
 !--            Upper left / front vertex
                IF ( point_index(k+1,j,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k+1,j,i) = point_index_count
                   point_index_count    = point_index_count + 1
@@ -781,6 +788,6 @@
                   surfaces%points(3,surfaces%npoints) = zw(k)
                ENDIF
-               
-               npg = npg + 1            
+
+               npg = npg + 1
                IF ( l == 0  .OR.  l == 1 )  THEN
                   surfaces%polygons(1,npg)   = 4
@@ -788,5 +795,5 @@
                   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)  
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
                ELSE
                   surfaces%polygons(1,npg)   = 4
@@ -796,15 +803,15 @@
                   surfaces%polygons(5,npg)   = point_index(k+1,j,i)
                ENDIF
-               
+
             ENDDO
-         
+
             DO  m = 1, surf_lsm_v(l)%ns
                i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
                j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
                k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
-!        
+!
 !--            Lower left /front vertex
                IF ( point_index(k,j,i) < 0 )  THEN
-                  surfaces%npoints   = surfaces%npoints + 1    
+                  surfaces%npoints   = surfaces%npoints + 1
                   point_index(k,j,i) = point_index_count
                   point_index_count  = point_index_count + 1
@@ -812,10 +819,10 @@
                   surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k-1)
-               ENDIF 
-!        
+               ENDIF
+!
 !--            Upper / lower right index for north- and south-facing surfaces
                IF ( l == 0  .OR.  l == 1 )  THEN
                   IF ( point_index(k,j,i+1) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j,i+1) = point_index_count
                      point_index_count    = point_index_count + 1
@@ -823,7 +830,7 @@
                      surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
                      surfaces%points(3,surfaces%npoints) = zw(k-1)
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j,i+1) = point_index_count
                      point_index_count      = point_index_count + 1
@@ -832,9 +839,9 @@
                      surfaces%points(3,surfaces%npoints) = zw(k)
                   ENDIF
-!        
+!
 !--            Upper / lower front index for east- and west-facing surfaces
                ELSEIF ( l == 2  .OR.  l == 3 )  THEN
                   IF ( point_index(k,j+1,i) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j+1,i) = point_index_count
                      point_index_count    = point_index_count + 1
@@ -842,7 +849,7 @@
                      surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
                      surfaces%points(3,surfaces%npoints) = zw(k-1)
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j+1,i) = point_index_count
                      point_index_count      = point_index_count + 1
@@ -852,8 +859,8 @@
                   ENDIF
                ENDIF
-!        
+!
 !--            Upper left / front vertex
                IF ( point_index(k+1,j,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k+1,j,i) = point_index_count
                   point_index_count    = point_index_count + 1
@@ -862,6 +869,6 @@
                   surfaces%points(3,surfaces%npoints) = zw(k)
                ENDIF
-               
-               npg = npg + 1            
+
+               npg = npg + 1
                IF ( l == 0  .OR.  l == 1 )  THEN
                   surfaces%polygons(1,npg)   = 4
@@ -869,5 +876,5 @@
                   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)  
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
                ELSE
                   surfaces%polygons(1,npg)   = 4
@@ -882,8 +889,8 @@
                j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
                k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
-!        
+!
 !--            Lower left /front vertex
                IF ( point_index(k,j,i) < 0 )  THEN
-                  surfaces%npoints   = surfaces%npoints + 1    
+                  surfaces%npoints   = surfaces%npoints + 1
                   point_index(k,j,i) = point_index_count
                   point_index_count  = point_index_count + 1
@@ -891,10 +898,10 @@
                   surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                   surfaces%points(3,surfaces%npoints) = zw(k-1)
-               ENDIF 
-!        
+               ENDIF
+!
 !--            Upper / lower right index for north- and south-facing surfaces
                IF ( l == 0  .OR.  l == 1 )  THEN
                   IF ( point_index(k,j,i+1) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j,i+1) = point_index_count
                      point_index_count    = point_index_count + 1
@@ -902,7 +909,7 @@
                      surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
                      surfaces%points(3,surfaces%npoints) = zw(k-1)
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j,i+1) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j,i+1) = point_index_count
                      point_index_count      = point_index_count + 1
@@ -911,9 +918,9 @@
                      surfaces%points(3,surfaces%npoints) = zw(k)
                   ENDIF
-!        
+!
 !--            Upper / lower front index for east- and west-facing surfaces
                ELSEIF ( l == 2  .OR.  l == 3 )  THEN
                   IF ( point_index(k,j+1,i) < 0 )  THEN
-                     surfaces%npoints     = surfaces%npoints + 1    
+                     surfaces%npoints     = surfaces%npoints + 1
                      point_index(k,j+1,i) = point_index_count
                      point_index_count    = point_index_count + 1
@@ -921,7 +928,7 @@
                      surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
                      surfaces%points(3,surfaces%npoints) = zw(k-1)
-                  ENDIF 
+                  ENDIF
                   IF ( point_index(k+1,j+1,i) < 0 )  THEN
-                     surfaces%npoints       = surfaces%npoints + 1    
+                     surfaces%npoints       = surfaces%npoints + 1
                      point_index(k+1,j+1,i) = point_index_count
                      point_index_count      = point_index_count + 1
@@ -931,8 +938,8 @@
                   ENDIF
                ENDIF
-!        
+!
 !--            Upper left / front vertex
                IF ( point_index(k+1,j,i) < 0 )  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
+                  surfaces%npoints     = surfaces%npoints + 1
                   point_index(k+1,j,i) = point_index_count
                   point_index_count    = point_index_count + 1
@@ -941,6 +948,6 @@
                   surfaces%points(3,surfaces%npoints) = zw(k)
                ENDIF
-               
-               npg = npg + 1            
+
+               npg = npg + 1
                IF ( l == 0  .OR.  l == 1 )  THEN
                   surfaces%polygons(1,npg)   = 4
@@ -948,5 +955,5 @@
                   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)  
+                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
                ELSE
                   surfaces%polygons(1,npg)   = 4
@@ -957,12 +964,12 @@
                ENDIF
             ENDDO
-         
+
          ENDDO
-!        
+!
 !--      Deallocate temporary dummy variable
          DEALLOCATE ( point_index )
 !
 !--      Sum-up total number of vertices on domain. This
-!--      will be needed for post-processing. 
+!--      will be needed for post-processing.
          surfaces%npoints_total = 0
 #if defined( __parallel )
@@ -974,7 +981,7 @@
        ENDIF
 !
-!--    If output to netcdf is enabled, set-up the coordinate arrays that 
-!--    unambiguously describe the position and orientation of each surface 
-!--    element. 
+!--    If output to netcdf is enabled, set-up the coordinate arrays that
+!--    unambiguously describe the position and orientation of each surface
+!--    element.
        IF ( to_netcdf )  THEN
 !
@@ -990,5 +997,5 @@
 !
 !--       Gather the number of surface on each processor, in order to number
-!--       the surface elements in ascending order with respect to the total 
+!--       the surface elements in ascending order with respect to the total
 !--       number of surfaces in the domain.
 #if defined( __parallel )
@@ -999,5 +1006,5 @@
 #endif
 !
-!--       First, however, determine the offset where couting of the surfaces 
+!--       First, however, determine the offset where couting of the surfaces
 !--       should start (the sum of surfaces on all PE's with lower MPI rank).
           i           = 0
@@ -1007,7 +1014,7 @@
              i           = i + 1
           ENDDO
-!       
-!--       Set coordinate arrays. For horizontal surfaces, azimuth 
-!--       angles are not defined (fill value). Zenith angle is 0 (180) for 
+!
+!--       Set coordinate arrays. For horizontal surfaces, azimuth
+!--       angles are not defined (fill value). Zenith angle is 0 (180) for
 !--       upward (downward)-facing surfaces.
           i  = start_count
@@ -1053,8 +1060,8 @@
              mm                   = mm + 1
           ENDDO
-!       
-!--       For vertical surfaces, zenith angles are not defined (fill value).  
-!--       Azimuth angle: eastward (0), westward (180), northward (270), 
-!--       southward (90). 
+!
+!--       For vertical surfaces, zenith angles are not defined (fill value).
+!--       Azimuth angle: eastward (0), westward (180), northward (270),
+!--       southward (90).
           DO  l = 0, 3
              IF ( l == 0 )  THEN
@@ -1075,5 +1082,5 @@
                 off_y =  0.0_wp
              ENDIF
-                
+
              DO  m = 1, surf_def_v(l)%ns
                 surfaces%s(mm)       = i
@@ -1107,12 +1114,12 @@
              ENDDO
           ENDDO
-!       
-!--       Finally, define UTM coordinates, which are the x/y-coordinates 
+!
+!--       Finally, define UTM coordinates, which are the x/y-coordinates
 !--       plus the origin (lower-left coordinate of the model domain).
-          surfaces%es_utm = surfaces%xs + init_model%origin_x 
-          surfaces%ns_utm = surfaces%ys + init_model%origin_y 
-!
-!--       Initialize NetCDF data output. Please note, local start position for 
-!--       the surface elements in the NetCDF file is surfaces%s(1), while 
+          surfaces%es_utm = surfaces%xs + init_model%origin_x
+          surfaces%ns_utm = surfaces%ys + init_model%origin_y
+!
+!--       Initialize NetCDF data output. Please note, local start position for
+!--       the surface elements in the NetCDF file is surfaces%s(1), while
 !--       the number of surfaces on the subdomain is given by surfaces%ns.
 #if defined( __netcdf4_parallel )
@@ -1120,15 +1127,17 @@
 !
 !--       Calculate number of time steps to be output
-          ntdim_surf(0) = dosurf_time_count(0) + CEILING(                            &
-                        ( end_time - MAX(                                            &
-                            MERGE(skip_time_dosurf, skip_time_dosurf + spinup_time,  &
-                                  data_output_during_spinup ),                       &
-                            simulated_time_at_begin )                                &
+          ntdim_surf(0) = dosurf_time_count(0) + CEILING(                      &
+                        ( end_time - MAX(                                      &
+                            MERGE( skip_time_dosurf,                           &
+                                   skip_time_dosurf + spinup_time,             &
+                                   data_output_during_spinup ),                &
+                            simulated_time_at_begin )                          &
                         ) / dt_dosurf )
 
           ntdim_surf(1) = dosurf_time_count(1) + CEILING(                      &
                         ( end_time - MAX(                                      &
-                            MERGE(   skip_time_dosurf_av, skip_time_dosurf_av  &
-                                  + spinup_time, data_output_during_spinup ),  &
+                            MERGE( skip_time_dosurf_av,                        &
+                                   skip_time_dosurf_av + spinup_time,          &
+                                   data_output_during_spinup ),                &
                             simulated_time_at_begin )                          &
                         ) / dt_dosurf_av )
@@ -1152,15 +1161,24 @@
                                     IOR( NF90_NOCLOBBER,                       &
                                     IOR( NF90_NETCDF4, NF90_MPIIO ) ),         &
-                                    id_set_surf(av), COMM = comm2d, INFO = MPI_INFO_NULL )
+                                    id_set_surf(av),                           &
+                                    COMM = comm2d, INFO = MPI_INFO_NULL )
              CALL netcdf_handle_error( 'surface_data_output_mod', 5550 )
 
              !- Write some global attributes
              IF ( av == 0 )  THEN
-                CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data', TRIM( run_description_header ), 5551 )
+                CALL netcdf_create_global_atts( id_set_surf(av),               &
+                                                'surface-data',                &
+                                                TRIM( run_description_header ),&
+                                                5551 )
                 time_average_text = ' '
              ELSE
-                CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data_av', TRIM( run_description_header ), 5552 )
-                WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_surf
-                nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'time_avg',  &
+                CALL netcdf_create_global_atts( id_set_surf(av),               &
+                                                'surface-data_av',             &
+                                                TRIM( run_description_header ),&
+                                                5552 )
+                WRITE ( time_average_text,'(F7.1,'' s avg'')' )  &
+                   averaging_interval_surf
+                nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL,          &
+                                        'time_avg',                            &
                                         TRIM( time_average_text ) )
                 CALL netcdf_handle_error( 'surface_data_output_mod', 5553 )
@@ -1170,15 +1188,22 @@
 !
 !--          Define time coordinate for surface data.
-!--          For parallel output the time dimension has to be limited (ntdim_surf),
-!--          otherwise the performance drops significantly.
+!--          For parallel output the time dimension has to be limited
+!--          (ntdim_surf), otherwise the performance drops significantly.
              CALL netcdf_create_dim( id_set_surf(av), 'time', ntdim_surf(av),  &
                                      id_dim_time_surf(av), 5554 )
 
-             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_time_surf(av) /),    &
-                                     'time', NF90_DOUBLE, id_var_time_surf(av),      &
-                                     'seconds since '//TRIM(init_model%origin_time), &
+             CALL netcdf_create_var( id_set_surf(av),                          &
+                                     (/ id_dim_time_surf(av) /),               &
+                                     'time', NF90_DOUBLE,                      &
+                                     id_var_time_surf(av),                     &
+                                     'seconds since '//                        &
+                                     TRIM(init_model%origin_time),             &
                                      'time', 5555, 5555, 5555 )
-             CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'standard_name', 'time', 5556)
-             CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'axis', 'T', 5557)
+
+             CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av),    &
+                                     'standard_name', 'time', 5556)
+
+             CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av),    &
+                                     'axis', 'T', 5557)
 !
 !--          Define spatial dimensions and coordinates:
@@ -1188,29 +1213,56 @@
              CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
                                      's', NF90_DOUBLE, id_var_s_surf(av),      &
-                                     '1', '', 5559, 5559, 5559 )
+                                     '1', 'number of surface element',         &
+                                     5559, 5559, 5559 )
 !
 !--          Define x coordinate
              CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
                                      'xs', NF90_DOUBLE, id_var_xs_surf(av),    &
-                                     'meters', '', 5561, 5561, 5561 )
+                                     'meters',                                 &
+                                     'distance to origin in x-direction',      &
+                                     5561, 5561, 5561 )
 !
 !--           Define y coordinate
              CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
                                      'ys', NF90_DOUBLE, id_var_ys_surf(av),    &
-                                     'meters', '', 5562, 5562, 5562 )
+                                     'meters',                                 &
+                                     'distance to origin in y-direction',      &
+                                     5562, 5562, 5562 )
 !
 !--          Define z coordinate
              CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
                                      'zs', NF90_DOUBLE, id_var_zs_surf(av),    &
-                                     'meters', '', 5560, 5560, 5560 )
+                                     'meters', 'height', 5560, 5560, 5560 )
+             CALL netcdf_create_att( id_set_surf(av), id_var_zs_surf(av),      &
+                                     'standard_name', 'height', 5583 )
 
 !
 !--          Define UTM coordinates
-             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /),    &
-                                     'Es_UTM', NF90_DOUBLE, id_var_etum_surf(av), &
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
+                                     'Es_UTM', NF90_DOUBLE,                    &
+                                     id_var_etum_surf(av),                     &
                                      'meters', '', 5563, 5563, 5563 )
-             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /),    &
-                                     'Ns_UTM', NF90_DOUBLE, id_var_nutm_surf(av), &
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
+                                     'Ns_UTM', NF90_DOUBLE,                    &
+                                     id_var_nutm_surf(av),                     &
                                      'meters', '', 5564, 5564, 5564 )
+
+!
+!--          Define angles
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
+                                     'azimuth', NF90_DOUBLE,                   &
+                                     id_var_azimuth_surf(av),                  &
+                                     'degree', 'azimuth angle',                &
+                                     5577, 5578, 5579,                         &
+                                     fill = .TRUE. )
+             CALL netcdf_create_att( id_set_surf(av), id_var_azimuth_surf(av), &
+                                     'standard_name', 'surface_azimuth_angle', &
+                                     5584 )
+
+             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
+                                     'zenith', NF90_DOUBLE,                    &
+                                     id_var_zenith_surf(av),                   &
+                                     'degree', '', 5580, 5581, 5582,           &
+                                     fill = .TRUE. )
 !
 !--          Define the variables
@@ -1227,12 +1279,12 @@
 !
 !--                Set no fill for every variable to increase performance.
-                nc_stat = NF90_DEF_VAR_FILL( id_set_surf(av),     &
-                                             id_var_dosurf(av,i), &
+                nc_stat = NF90_DEF_VAR_FILL( id_set_surf(av),                  &
+                                             id_var_dosurf(av,i),              &
                                              1, 0 )
                 CALL netcdf_handle_error( 'surface_data_output_init', 5566 )
 !
 !--                Set collective io operations for parallel io
-                nc_stat = NF90_VAR_PAR_ACCESS( id_set_surf(av),     &
-                                               id_var_dosurf(av,i), &
+                nc_stat = NF90_VAR_PAR_ACCESS( id_set_surf(av),                &
+                                               id_var_dosurf(av,i),            &
                                                NF90_COLLECTIVE )
                 CALL netcdf_handle_error( 'surface_data_output_init', 5567 )
@@ -1250,6 +1302,6 @@
 
 !
-!--          Set general no fill, otherwise the performance drops significantly for
-!--          parallel output.
+!--          Set general no fill, otherwise the performance drops significantly
+!--          for parallel output.
              nc_stat = NF90_SET_FILL( id_set_surf(av), NF90_NOFILL, oldmode )
              CALL netcdf_handle_error( 'surface_data_output_init', 5569 )
@@ -1272,6 +1324,6 @@
                 ENDDO
 
-                nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_s_surf(av),  &
-                                        netcdf_data_1d, start = (/ 1 /),     &
+                nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_s_surf(av),    &
+                                        netcdf_data_1d, start = (/ 1 /),       &
                                         count = (/ surfaces%ns_total /) )
                 CALL netcdf_handle_error( 'surface_data_output_init', 5571 )
@@ -1298,14 +1350,28 @@
              CALL netcdf_handle_error( 'surface_data_output_init', 5574 )
 
-             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_etum_surf(av),        &
-                                     surfaces%es_utm, start = (/ surfaces%s(1) /), &
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_etum_surf(av),    &
+                                     surfaces%es_utm,                          &
+                                     start = (/ surfaces%s(1) /),              &
                                      count = (/ surfaces%ns /) )
              CALL netcdf_handle_error( 'surface_data_output_init', 5575 )
 
-             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_nutm_surf(av),        &
-                                     surfaces%ns_utm, start = (/ surfaces%s(1) /), &
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_nutm_surf(av),    &
+                                     surfaces%ns_utm,                          &
+                                     start = (/ surfaces%s(1) /),              &
                                      count = (/ surfaces%ns /) )
              CALL netcdf_handle_error( 'surface_data_output_init', 5576 )
-             
+
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_azimuth_surf(av), &
+                                     surfaces%azimuth,                         &
+                                     start = (/ surfaces%s(1) /),              &
+                                     count = (/ surfaces%ns /) )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5585 )
+
+             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_zenith_surf(av),  &
+                                     surfaces%zenith,                          &
+                                     start = (/ surfaces%s(1) /),              &
+                                     count = (/ surfaces%ns /) )
+             CALL netcdf_handle_error( 'surface_data_output_init', 5586 )
+
           ENDDO
 #endif
@@ -1314,15 +1380,15 @@
 
    END SUBROUTINE surface_data_output_init
-   
+
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Routine for controlling the data output. Surface data is collected from 
-!> different types of surfaces (default, natural, urban) and different 
+!> 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_data_output( av )
-   
+
       USE control_parameters,                                                  &
           ONLY:  io_blocks, io_group, time_since_reference_point
@@ -1331,9 +1397,9 @@
           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
@@ -1344,5 +1410,5 @@
       IF ( dosurf_no(av) == 0 )  RETURN
 !
-!--   In case of VTK output, check if binary files are open and write coordinates. 
+!--   In case of VTK output, check if binary files are open and write coordinates.
       IF ( to_vtk )  THEN
 
@@ -1372,9 +1438,10 @@
 #if defined( __netcdf4_parallel )
          IF ( dosurf_time_count(av) + 1 > ntdim_surf(av) )  THEN
-            WRITE ( message_string, * ) 'Output of surface data is not given at t=',          &
-                                        time_since_reference_point, 's because the maximum ', & 
-                                        'number of output time levels is exceeded.'
+            WRITE ( message_string, * )                               &
+               'Output of surface data is not given at t=',           &
+               time_since_reference_point, 's because the maximum ',  &
+               'number of output time levels is exceeded.'
             CALL message( 'surface_data_output', 'PA0539', 0, 1, 0, 6, 0 )
-            
+
             RETURN
 
@@ -1403,5 +1470,5 @@
          trimvar = TRIM( dosurf(av,n_out) )
 !
-!--      Set the output array to the _FillValue in case it is not 
+!--      Set the output array to the _FillValue in case it is not
 !--      defined for each type of surface.
          surfaces%var_out = surfaces%fillvalue
@@ -1427,12 +1494,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1456,12 +1523,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1485,12 +1552,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1514,12 +1581,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1543,12 +1610,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1572,12 +1639,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1601,12 +1668,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1630,12 +1697,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1659,12 +1726,12 @@
                                                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
-                                                        
+                                               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
 
@@ -1695,5 +1762,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -1724,5 +1791,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -1746,5 +1813,5 @@
                                                surf_def_v(3)%z0q,              &
                                                surf_lsm_v(3)%z0q,              &
-                                               surf_usm_v(3)%z0q ) 
+                                               surf_usm_v(3)%z0q )
                ELSE
 !
@@ -1782,5 +1849,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -1811,5 +1878,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -1840,5 +1907,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -1869,5 +1936,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -1898,5 +1965,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -1955,5 +2022,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -1984,5 +2051,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2013,5 +2080,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2042,5 +2109,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2071,5 +2138,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2100,5 +2167,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2129,5 +2196,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2158,5 +2225,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2187,5 +2254,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2216,5 +2283,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2238,12 +2305,12 @@
                                                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
-                                                        
+                                               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
 
@@ -2274,5 +2341,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2303,5 +2370,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2332,5 +2399,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2361,5 +2428,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2383,5 +2450,5 @@
                                                surf_def_v(3)%ghf,              &
                                                surf_lsm_v(3)%ghf,              &
-                                               surf_usm_v(3)%ghf )      
+                                               surf_usm_v(3)%ghf )
                                                                         ELSE
 !
@@ -2390,8 +2457,8 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
-               ENDIF
-                      
-            CASE ( 'r_a' )                                                     
+
+               ENDIF
+
+            CASE ( 'r_a' )
 !
 !--            Output of instantaneous data
@@ -2412,5 +2479,5 @@
                                                surf_def_v(3)%r_a,              &
                                                surf_lsm_v(3)%r_a,              &
-                                               surf_usm_v(3)%r_a )      
+                                               surf_usm_v(3)%r_a )
                                                                       ELSE
 !
@@ -2419,8 +2486,8 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
-               ENDIF
-                        
-            CASE ( 'r_soil' )                                                  
+
+               ENDIF
+
+            CASE ( 'r_soil' )
 !
 !--            Output of instantaneous data
@@ -2441,12 +2508,12 @@
                                                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
-                                                        
+                                               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
 
@@ -2470,12 +2537,12 @@
                                                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
-                                                        
+                                               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
 
@@ -2499,12 +2566,12 @@
                                                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
-                                                        
+                                               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
 
@@ -2535,5 +2602,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2564,5 +2631,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2593,5 +2660,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2622,5 +2689,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2651,5 +2718,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2680,5 +2747,5 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
+
                ENDIF
 
@@ -2709,7 +2776,7 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
-               ENDIF
-               
+
+               ENDIF
+
             CASE ( 'uvw1' )
 !
@@ -2738,6 +2805,6 @@
                                         REAL( average_count_surf, KIND=wp )
                   surfaces%var_av(:,n_out) = 0.0_wp
-                                                        
-               ENDIF   
+
+               ENDIF
 
 !
@@ -2747,5 +2814,5 @@
          END SELECT
 !
-!--      Write to binary file: 
+!--      Write to binary file:
 !--      - surfaces%points ( 3, 1-npoints )
 !--      - surfaces%polygons ( 5, 1-ns )
@@ -2768,12 +2835,13 @@
             ENDDO
          ENDIF
-         
+
          IF ( to_netcdf )  THEN
 #if defined( __netcdf4_parallel )
 !
 !--         Write output array to file
-            nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_dosurf(av,n_out),            &
-                                    surfaces%var_out,                                    &
-                                    start = (/ surfaces%s(1), dosurf_time_count(av) /),  &
+            nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_dosurf(av,n_out),  &
+                                    surfaces%var_out,                          &
+                                    start = (/ surfaces%s(1),                  &
+                                               dosurf_time_count(av) /),       &
                                     count = (/ surfaces%ns, 1 /) )
             CALL netcdf_handle_error( 'surface_data_output', 6667 )
@@ -2782,27 +2850,27 @@
 
       ENDDO
-      
+
 !
 !--   If averaged output was written to NetCDF file, set the counter to zero
       IF ( av == 1 )  average_count_surf = 0
-             
+
    END SUBROUTINE surface_data_output
-   
+
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Routine for controlling the data averaging. 
-!------------------------------------------------------------------------------!     
+!> Routine for controlling the data averaging.
+!------------------------------------------------------------------------------!
    SUBROUTINE surface_data_output_averaging
-   
+
       IMPLICIT NONE
-      
+
       CHARACTER(LEN=100) ::  trimvar !< dummy variable for current output variable
-       
+
       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) )
@@ -2826,5 +2894,5 @@
                                            surf_def_v(3)%us,                   &
                                            surf_lsm_v(3)%us,                   &
-                                           surf_usm_v(3)%us, n_out ) 
+                                           surf_usm_v(3)%us, n_out )
 
             CASE ( 'ts' )
@@ -2844,5 +2912,5 @@
                                            surf_def_v(3)%ts,                   &
                                            surf_lsm_v(3)%ts,                   &
-                                           surf_usm_v(3)%ts, n_out ) 
+                                           surf_usm_v(3)%ts, n_out )
 
             CASE ( 'qs' )
@@ -2862,5 +2930,5 @@
                                            surf_def_v(3)%qs,                   &
                                            surf_lsm_v(3)%qs,                   &
-                                           surf_usm_v(3)%qs, n_out ) 
+                                           surf_usm_v(3)%qs, n_out )
 
             CASE ( 'ss' )
@@ -2880,5 +2948,5 @@
                                            surf_def_v(3)%ss,                   &
                                            surf_lsm_v(3)%ss,                   &
-                                           surf_usm_v(3)%ss, n_out ) 
+                                           surf_usm_v(3)%ss, n_out )
 
             CASE ( 'qcs' )
@@ -2898,5 +2966,5 @@
                                            surf_def_v(3)%qcs,                  &
                                            surf_lsm_v(3)%qcs,                  &
-                                           surf_usm_v(3)%qcs, n_out ) 
+                                           surf_usm_v(3)%qcs, n_out )
 
             CASE ( 'ncs' )
@@ -2916,5 +2984,5 @@
                                            surf_def_v(3)%ncs,                  &
                                            surf_lsm_v(3)%ncs,                  &
-                                           surf_usm_v(3)%ncs, n_out ) 
+                                           surf_usm_v(3)%ncs, n_out )
 
             CASE ( 'qrs' )
@@ -2934,5 +3002,5 @@
                                            surf_def_v(3)%qrs,                  &
                                            surf_lsm_v(3)%qrs,                  &
-                                           surf_usm_v(3)%qrs, n_out ) 
+                                           surf_usm_v(3)%qrs, n_out )
 
             CASE ( 'nrs' )
@@ -2952,5 +3020,5 @@
                                            surf_def_v(3)%nrs,                  &
                                            surf_lsm_v(3)%nrs,                  &
-                                           surf_usm_v(3)%nrs, n_out ) 
+                                           surf_usm_v(3)%nrs, n_out )
 
             CASE ( 'ol' )
@@ -2970,5 +3038,5 @@
                                            surf_def_v(3)%ol,                   &
                                            surf_lsm_v(3)%ol,                   &
-                                           surf_usm_v(3)%ol, n_out ) 
+                                           surf_usm_v(3)%ol, n_out )
 
             CASE ( 'z0' )
@@ -2988,5 +3056,5 @@
                                            surf_def_v(3)%z0,                   &
                                            surf_lsm_v(3)%z0,                   &
-                                           surf_usm_v(3)%z0, n_out ) 
+                                           surf_usm_v(3)%z0, n_out )
 
             CASE ( 'z0h' )
@@ -3006,5 +3074,5 @@
                                            surf_def_v(3)%z0h,                  &
                                            surf_lsm_v(3)%z0h,                  &
-                                           surf_usm_v(3)%z0h, n_out )          
+                                           surf_usm_v(3)%z0h, n_out )
 
             CASE ( 'z0q' )
@@ -3024,7 +3092,7 @@
                                            surf_def_v(3)%z0q,                  &
                                            surf_lsm_v(3)%z0q,                  &
-                                           surf_usm_v(3)%z0q, n_out ) 
-
-            CASE ( 'theta1' )                                                  
+                                           surf_usm_v(3)%z0q, n_out )
+
+            CASE ( 'theta1' )
                CALL surface_data_output_sum_up( surf_def_h(0)%pt1,             &
                                            surf_def_h(1)%pt1,                  &
@@ -3042,7 +3110,7 @@
                                            surf_def_v(3)%pt1,                  &
                                            surf_lsm_v(3)%pt1,                  &
-                                           surf_usm_v(3)%pt1, n_out )          
-                                                                               
-            CASE ( 'qv1' )                                                     
+                                           surf_usm_v(3)%pt1, n_out )
+
+            CASE ( 'qv1' )
                CALL surface_data_output_sum_up( surf_def_h(0)%qv1,             &
                                            surf_def_h(1)%qv1,                  &
@@ -3060,5 +3128,5 @@
                                            surf_def_v(3)%qv1,                  &
                                            surf_lsm_v(3)%qv1,                  &
-                                           surf_usm_v(3)%qv1, n_out ) 
+                                           surf_usm_v(3)%qv1, n_out )
 
             CASE ( 'thetav1' )
@@ -3078,5 +3146,5 @@
                                            surf_def_v(3)%vpt1,                 &
                                            surf_lsm_v(3)%vpt1,                 &
-                                           surf_usm_v(3)%vpt1, n_out ) 
+                                           surf_usm_v(3)%vpt1, n_out )
 
             CASE ( 'usws' )
@@ -3096,5 +3164,5 @@
                                            surf_def_v(3)%usws,                 &
                                            surf_lsm_v(3)%usws,                 &
-                                           surf_usm_v(3)%usws, n_out ) 
+                                           surf_usm_v(3)%usws, n_out )
 
             CASE ( 'vsws' )
@@ -3114,5 +3182,5 @@
                                            surf_def_v(3)%vsws,                 &
                                            surf_lsm_v(3)%vsws,                 &
-                                           surf_usm_v(3)%vsws, n_out ) 
+                                           surf_usm_v(3)%vsws, n_out )
 
             CASE ( 'shf' )
@@ -3120,5 +3188,5 @@
                                            surf_def_h(1)%shf,                  &
                                            surf_lsm_h%shf,                     &
-                                           surf_usm_h%shf,                     &  
+                                           surf_usm_h%shf,                     &
                                            surf_def_v(0)%shf,                  &
                                            surf_lsm_v(0)%shf,                  &
@@ -3150,5 +3218,5 @@
                                            surf_def_v(3)%qsws,                 &
                                            surf_lsm_v(3)%qsws,                 &
-                                           surf_usm_v(3)%qsws, n_out ) 
+                                           surf_usm_v(3)%qsws, n_out )
 
             CASE ( 'ssws' )
@@ -3168,5 +3236,5 @@
                                            surf_def_v(3)%ssws,                 &
                                            surf_lsm_v(3)%ssws,                 &
-                                           surf_usm_v(3)%ssws, n_out ) 
+                                           surf_usm_v(3)%ssws, n_out )
 
             CASE ( 'qcsws' )
@@ -3186,5 +3254,5 @@
                                            surf_def_v(3)%qcsws,                &
                                            surf_lsm_v(3)%qcsws,                &
-                                           surf_usm_v(3)%qcsws, n_out ) 
+                                           surf_usm_v(3)%qcsws, n_out )
 
             CASE ( 'ncsws' )
@@ -3204,5 +3272,5 @@
                                            surf_def_v(3)%ncsws,                &
                                            surf_lsm_v(3)%ncsws,                &
-                                           surf_usm_v(3)%ncsws, n_out ) 
+                                           surf_usm_v(3)%ncsws, n_out )
 
             CASE ( 'qrsws' )
@@ -3222,5 +3290,5 @@
                                            surf_def_v(3)%qrsws,                &
                                            surf_lsm_v(3)%qrsws,                &
-                                           surf_usm_v(3)%qrsws, n_out ) 
+                                           surf_usm_v(3)%qrsws, n_out )
 
             CASE ( 'nrsws' )
@@ -3240,5 +3308,5 @@
                                            surf_def_v(3)%nrsws,                &
                                            surf_lsm_v(3)%nrsws,                &
-                                           surf_usm_v(3)%nrsws, n_out ) 
+                                           surf_usm_v(3)%nrsws, n_out )
 
             CASE ( 'sasws' )
@@ -3258,5 +3326,5 @@
                                            surf_def_v(3)%sasws,                &
                                            surf_lsm_v(3)%sasws,                &
-                                           surf_usm_v(3)%sasws, n_out ) 
+                                           surf_usm_v(3)%sasws, n_out )
 
             CASE ( 'q_surface' )
@@ -3276,6 +3344,6 @@
                                            surf_def_v(3)%q_surface,            &
                                            surf_lsm_v(3)%q_surface,            &
-                                           surf_usm_v(3)%q_surface, n_out ) 
-                                           
+                                           surf_usm_v(3)%q_surface, n_out )
+
 
             CASE ( 'theta_surface' )
@@ -3295,5 +3363,5 @@
                                            surf_def_v(3)%pt_surface,           &
                                            surf_lsm_v(3)%pt_surface,           &
-                                           surf_usm_v(3)%pt_surface, n_out ) 
+                                           surf_usm_v(3)%pt_surface, n_out )
 
             CASE ( 'thetav_surface' )
@@ -3313,5 +3381,5 @@
                                            surf_def_v(3)%vpt_surface,          &
                                            surf_lsm_v(3)%vpt_surface,          &
-                                           surf_usm_v(3)%vpt_surface, n_out ) 
+                                           surf_usm_v(3)%vpt_surface, n_out )
 
             CASE ( 'rad_net' )
@@ -3331,5 +3399,5 @@
                                            surf_def_v(3)%rad_net,              &
                                            surf_lsm_v(3)%rad_net,              &
-                                           surf_usm_v(3)%rad_net, n_out ) 
+                                           surf_usm_v(3)%rad_net, n_out )
 
             CASE ( 'rad_lw_in' )
@@ -3349,5 +3417,5 @@
                                            surf_def_v(3)%rad_lw_in,            &
                                            surf_lsm_v(3)%rad_lw_in,            &
-                                           surf_usm_v(3)%rad_lw_in, n_out ) 
+                                           surf_usm_v(3)%rad_lw_in, n_out )
 
             CASE ( 'rad_lw_out' )
@@ -3367,5 +3435,5 @@
                                            surf_def_v(3)%rad_lw_out,           &
                                            surf_lsm_v(3)%rad_lw_out,           &
-                                           surf_usm_v(3)%rad_lw_out, n_out ) 
+                                           surf_usm_v(3)%rad_lw_out, n_out )
 
             CASE ( 'rad_sw_in' )
@@ -3385,5 +3453,5 @@
                                            surf_def_v(3)%rad_sw_in,            &
                                            surf_lsm_v(3)%rad_sw_in,            &
-                                           surf_usm_v(3)%rad_sw_in, n_out ) 
+                                           surf_usm_v(3)%rad_sw_in, n_out )
 
             CASE ( 'rad_sw_out' )
@@ -3403,5 +3471,5 @@
                                            surf_def_v(3)%rad_sw_out,           &
                                            surf_lsm_v(3)%rad_sw_out,           &
-                                           surf_usm_v(3)%rad_sw_out, n_out ) 
+                                           surf_usm_v(3)%rad_sw_out, n_out )
 
             CASE ( 'ghf' )
@@ -3421,7 +3489,7 @@
                                            surf_def_v(3)%ghf,                  &
                                            surf_lsm_v(3)%ghf,                  &
-                                           surf_usm_v(3)%ghf, n_out )          
-                                                                               
-            CASE ( 'r_a' )                                                     
+                                           surf_usm_v(3)%ghf, n_out )
+
+            CASE ( 'r_a' )
                CALL surface_data_output_sum_up( surf_def_h(0)%r_a,             &
                                            surf_def_h(1)%r_a,                  &
@@ -3439,7 +3507,7 @@
                                            surf_def_v(3)%r_a,                  &
                                            surf_lsm_v(3)%r_a,                  &
-                                           surf_usm_v(3)%r_a, n_out )          
-                                                                               
-            CASE ( 'r_soil' )                                                  
+                                           surf_usm_v(3)%r_a, n_out )
+
+            CASE ( 'r_soil' )
                CALL surface_data_output_sum_up( surf_def_h(0)%r_soil,          &
                                            surf_def_h(1)%r_soil,               &
@@ -3457,5 +3525,5 @@
                                            surf_def_v(3)%r_soil,               &
                                            surf_lsm_v(3)%r_soil,               &
-                                           surf_usm_v(3)%r_soil, n_out ) 
+                                           surf_usm_v(3)%r_soil, n_out )
 
             CASE ( 'r_canopy' )
@@ -3475,5 +3543,5 @@
                                            surf_def_v(3)%r_canopy,             &
                                            surf_lsm_v(3)%r_canopy,             &
-                                           surf_usm_v(3)%r_canopy, n_out ) 
+                                           surf_usm_v(3)%r_canopy, n_out )
 
             CASE ( 'r_s' )
@@ -3493,5 +3561,5 @@
                                            surf_def_v(3)%r_s,                  &
                                            surf_lsm_v(3)%r_s,                  &
-                                           surf_usm_v(3)%r_s, n_out ) 
+                                           surf_usm_v(3)%r_s, n_out )
 
 
@@ -3548,5 +3616,5 @@
                                            surf_lsm_v(3)%rad_sw_ref,           &
                                            surf_usm_v(3)%rad_sw_ref, n_out )
-               
+
             CASE ( 'rad_sw_res' )
                CALL surface_data_output_sum_up( surf_def_h(0)%rad_sw_res,      &
@@ -3620,5 +3688,5 @@
                                            surf_lsm_v(3)%rad_lw_res,           &
                                            surf_usm_v(3)%rad_lw_res, n_out )
-                                           
+
             CASE ( 'uvw1' )
                CALL surface_data_output_sum_up( surf_def_h(0)%uvw_abs,         &
@@ -3637,17 +3705,17 @@
                                            surf_def_v(3)%uvw_abs,              &
                                            surf_lsm_v(3)%uvw_abs,              &
-                                           surf_usm_v(3)%uvw_abs, n_out )                               
+                                           surf_usm_v(3)%uvw_abs, n_out )
 
          END SELECT
       ENDDO
-      
-   
+
+
    END SUBROUTINE surface_data_output_averaging
-   
+
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Sum-up the surface data for average output variables.  
-!------------------------------------------------------------------------------!   
+!> Sum-up the surface data for average output variables.
+!------------------------------------------------------------------------------!
    SUBROUTINE surface_data_output_sum_up( var_def_h0, var_def_h1,              &
                                      var_lsm_h,  var_usm_h,                    &
@@ -3656,5 +3724,5 @@
                                      var_def_v2, var_lsm_v2, var_usm_v2,       &
                                      var_def_v3, var_lsm_v3, var_usm_v3, n_out )
-   
+
       IMPLICIT NONE
 
@@ -3662,5 +3730,5 @@
       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
@@ -3679,21 +3747,21 @@
       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. 
+
+!
+!--   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. 
+!--   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 
+!--   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. 
+!--   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 
+         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)      &
@@ -3701,8 +3769,8 @@
          ENDDO
       ELSE
-         n_surf = n_surf + surf_def_h(0)%ns  
+         n_surf = n_surf + surf_def_h(0)%ns
       ENDIF
       IF ( ALLOCATED( var_def_h1 ) )  THEN
-         DO  m = 1, surf_def_h(1)%ns 
+         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)      &
@@ -3713,5 +3781,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_h ) )  THEN
-         DO  m = 1, surf_lsm_h%ns 
+         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)      &
@@ -3722,5 +3790,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_h ) )  THEN
-         DO  m = 1, surf_usm_h%ns 
+         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)      &
@@ -3733,5 +3801,5 @@
 !--   Write northward-facing
       IF ( ALLOCATED( var_def_v0 ) )  THEN
-         DO  m = 1, surf_def_v(0)%ns 
+         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)      &
@@ -3742,5 +3810,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_v0 ) )  THEN
-         DO  m = 1, surf_lsm_v(0)%ns 
+         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)      &
@@ -3751,5 +3819,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_v0 ) )  THEN
-         DO  m = 1, surf_usm_v(0)%ns 
+         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)      &
@@ -3762,5 +3830,5 @@
 !--   Write southward-facing
       IF ( ALLOCATED( var_def_v1 ) )  THEN
-         DO  m = 1, surf_def_v(1)%ns 
+         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)      &
@@ -3771,5 +3839,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_v1 ) )  THEN
-         DO  m = 1, surf_lsm_v(1)%ns 
+         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)      &
@@ -3780,5 +3848,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_v1 ) )  THEN
-         DO  m = 1, surf_usm_v(1)%ns 
+         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)      &
@@ -3791,5 +3859,5 @@
 !--   Write eastward-facing
       IF ( ALLOCATED( var_def_v2 ) )  THEN
-         DO  m = 1, surf_def_v(2)%ns 
+         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)      &
@@ -3800,5 +3868,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_v2 ) )  THEN
-         DO  m = 1, surf_lsm_v(2)%ns 
+         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)      &
@@ -3809,5 +3877,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_v2 ) )  THEN
-         DO  m = 1, surf_usm_v(2)%ns 
+         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)      &
@@ -3820,5 +3888,5 @@
 !--   Write westward-facing
       IF ( ALLOCATED( var_def_v3 ) )  THEN
-         DO  m = 1, surf_def_v(3)%ns 
+         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)      &
@@ -3829,5 +3897,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_v3 ) )  THEN
-         DO  m = 1, surf_lsm_v(3)%ns 
+         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)      &
@@ -3838,5 +3906,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_v3 ) )  THEN
-         DO  m = 1, surf_usm_v(3)%ns 
+         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)      &
@@ -3846,7 +3914,7 @@
          n_surf = n_surf + surf_usm_v(3)%ns
       ENDIF
-   
+
    END SUBROUTINE surface_data_output_sum_up
-   
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -3859,11 +3927,11 @@
                                       var_def_v1, var_lsm_v1, var_usm_v1,      &
                                       var_def_v2, var_lsm_v2, var_usm_v2,      &
-                                      var_def_v3, var_lsm_v3, var_usm_v3 )                                               
-   
+                                      var_def_v3, var_lsm_v3, 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
@@ -3882,29 +3950,29 @@
       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. 
+
+!
+!--   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. 
+!--   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 
+!--   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. 
+!--   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 
+         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  
+         n_surf = n_surf + surf_def_h(0)%ns
       ENDIF
       IF ( ALLOCATED( var_def_h1 ) )  THEN
-         DO  m = 1, surf_def_h(1)%ns 
+         DO  m = 1, surf_def_h(1)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_def_h1(m)
@@ -3914,5 +3982,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_h ) )  THEN
-         DO  m = 1, surf_lsm_h%ns 
+         DO  m = 1, surf_lsm_h%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_lsm_h(m)
@@ -3922,5 +3990,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_h ) )  THEN
-         DO  m = 1, surf_usm_h%ns 
+         DO  m = 1, surf_usm_h%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_usm_h(m)
@@ -3932,5 +4000,5 @@
 !--   Write northward-facing
       IF ( ALLOCATED( var_def_v0 ) )  THEN
-         DO  m = 1, surf_def_v(0)%ns 
+         DO  m = 1, surf_def_v(0)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_def_v0(m)
@@ -3940,5 +4008,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_v0 ) )  THEN
-         DO  m = 1, surf_lsm_v(0)%ns 
+         DO  m = 1, surf_lsm_v(0)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_lsm_v0(m)
@@ -3948,5 +4016,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_v0 ) )  THEN
-         DO  m = 1, surf_usm_v(0)%ns 
+         DO  m = 1, surf_usm_v(0)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_usm_v0(m)
@@ -3958,5 +4026,5 @@
 !--   Write southward-facing
       IF ( ALLOCATED( var_def_v1 ) )  THEN
-         DO  m = 1, surf_def_v(1)%ns 
+         DO  m = 1, surf_def_v(1)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_def_v1(m)
@@ -3966,5 +4034,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_v1 ) )  THEN
-         DO  m = 1, surf_lsm_v(1)%ns 
+         DO  m = 1, surf_lsm_v(1)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_lsm_v1(m)
@@ -3974,5 +4042,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_v1 ) )  THEN
-         DO  m = 1, surf_usm_v(1)%ns 
+         DO  m = 1, surf_usm_v(1)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_usm_v1(m)
@@ -3984,5 +4052,5 @@
 !--   Write eastward-facing
       IF ( ALLOCATED( var_def_v2 ) )  THEN
-         DO  m = 1, surf_def_v(2)%ns 
+         DO  m = 1, surf_def_v(2)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_def_v2(m)
@@ -3992,5 +4060,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_v2 ) )  THEN
-         DO  m = 1, surf_lsm_v(2)%ns 
+         DO  m = 1, surf_lsm_v(2)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_lsm_v2(m)
@@ -4000,5 +4068,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_v2 ) )  THEN
-         DO  m = 1, surf_usm_v(2)%ns 
+         DO  m = 1, surf_usm_v(2)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_usm_v2(m)
@@ -4010,5 +4078,5 @@
 !--   Write westward-facing
       IF ( ALLOCATED( var_def_v3 ) )  THEN
-         DO  m = 1, surf_def_v(3)%ns 
+         DO  m = 1, surf_def_v(3)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_def_v3(m)
@@ -4018,5 +4086,5 @@
       ENDIF
       IF ( ALLOCATED( var_lsm_v3 ) )  THEN
-         DO  m = 1, surf_lsm_v(3)%ns 
+         DO  m = 1, surf_lsm_v(3)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_lsm_v3(m)
@@ -4026,5 +4094,5 @@
       ENDIF
       IF ( ALLOCATED( var_usm_v3 ) )  THEN
-         DO  m = 1, surf_usm_v(3)%ns 
+         DO  m = 1, surf_usm_v(3)%ns
             n_surf                   = n_surf + 1
             surfaces%var_out(n_surf) = var_usm_v3(m)
@@ -4033,7 +4101,7 @@
          n_surf = n_surf + surf_usm_v(3)%ns
       ENDIF
-   
+
    END SUBROUTINE surface_data_output_collect
-   
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4045,6 +4113,6 @@
        IMPLICIT NONE
 
-       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
-       
+       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
+
 
        NAMELIST /surface_data_output_parameters/                               &
@@ -4053,7 +4121,7 @@
                                   skip_time_dosurf, skip_time_dosurf_av,       &
                                   to_netcdf, to_vtk
-                                  
+
        line = ' '
-  
+
 !
 !--    Try to find the namelist
@@ -4075,18 +4143,18 @@
  10    BACKSPACE( 11 )
        READ( 11 , '(A)') line
-       CALL parin_fail_message( 'surface_data_output_parameters', line )      
-       
+       CALL parin_fail_message( 'surface_data_output_parameters', line )
+
  14    CONTINUE
-       
+
 
     END SUBROUTINE surface_data_output_parin
-    
-    
+
+
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Check the input parameters for consistency. Further pre-process the given 
+!> 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. 
+!> variables and map them onto internal output array.
 !------------------------------------------------------------------------------!
     SUBROUTINE surface_data_output_check_parameters
@@ -4095,5 +4163,5 @@
            ONLY:  averaging_interval, dt_data_output, initializing_actions,    &
                   message_string
-                  
+
        USE pegrid,                                                             &
            ONLY:  numprocs_previous_run
@@ -4125,5 +4193,5 @@
                         'PA0536', 1, 2, 0, 6, 0 )
        ENDIF
-       
+
 #if ! defined( __netcdf4_parallel )
 !
@@ -4147,5 +4215,5 @@
       ENDIF
 !
-!--   Count number of output variables and separate output strings for 
+!--   Count number of output variables and separate output strings for
 !--   average and non-average output variables.
       n_out = 0
@@ -4194,5 +4262,5 @@
             CASE ( 'usws', 'vsws' )
                unit = 'm2/s2'
-            
+
             CASE ( 'qcsws', 'ncsws', 'qrsws', 'nrsws', 'sasws' )
 
@@ -4238,6 +4306,6 @@
 
     END SUBROUTINE surface_data_output_check_parameters
-    
-    
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4246,5 +4314,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE surface_data_output_last_action( av )
-  
+
       USE control_parameters,                                                  &
           ONLY:  io_blocks, io_group
@@ -4254,5 +4322,5 @@
 
       IMPLICIT NONE
-      
+
       INTEGER(iwp) ::  av     !< id indicating average or non-average data output
       INTEGER(iwp) ::  i      !< loop index
@@ -4276,9 +4344,9 @@
             CALL MPI_BARRIER( comm2d, ierr )
 #endif
-         ENDDO    
+         ENDDO
       ENDIF
 
-    END SUBROUTINE surface_data_output_last_action    
-    
+    END SUBROUTINE surface_data_output_last_action
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4291,25 +4359,25 @@
        USE control_parameters,                                                 &
            ONLY: length, restart_string
-       
+
        IMPLICIT NONE
-       
+
        LOGICAL, INTENT(OUT)  ::  found !< flag indicating if variable was found
-       
+
        found = .TRUE.
-       
+
        SELECT CASE ( restart_string(1:length) )
-       
+
           CASE ( 'average_count_surf' )
              READ ( 13 )  average_count_surf
-       
+
           CASE DEFAULT
-       
+
              found = .FALSE.
-       
+
        END SELECT
 
 
     END SUBROUTINE surface_data_output_rrd_global
-    
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4324,5 +4392,5 @@
        USE control_parameters,                                                 &
            ONLY: length, restart_string
-           
+
        IMPLICIT NONE
 
@@ -4332,13 +4400,13 @@
        INTEGER(iwp)       ::  ns_h_on_file_usm  !< number of horizontal surface elements (urban type) on file
        INTEGER(iwp)       ::  nxlc              !< index of left boundary on current subdomain
-       INTEGER(iwp)       ::  nxlf              !< index of left boundary on former subdomain 
-       INTEGER(iwp)       ::  nxl_on_file       !< index of left boundary on former local domain 
+       INTEGER(iwp)       ::  nxlf              !< index of left boundary on former subdomain
+       INTEGER(iwp)       ::  nxl_on_file       !< index of left boundary on former local domain
        INTEGER(iwp)       ::  nxrc              !< index of right boundary on current subdomain
        INTEGER(iwp)       ::  nxrf              !< index of right boundary on former subdomain
-       INTEGER(iwp)       ::  nxr_on_file       !< index of right boundary on former local domain 
+       INTEGER(iwp)       ::  nxr_on_file       !< index of right boundary on former local domain
        INTEGER(iwp)       ::  nync              !< index of north boundary on current subdomain
        INTEGER(iwp)       ::  nynf              !< index of north boundary on former subdomain
-       INTEGER(iwp)       ::  nyn_on_file       !< index of north boundary on former local domain 
-       INTEGER(iwp)       ::  nysc              !< index of south boundary on current subdomain 
+       INTEGER(iwp)       ::  nyn_on_file       !< index of north boundary on former local domain
+       INTEGER(iwp)       ::  nysc              !< index of south boundary on current subdomain
        INTEGER(iwp)       ::  nysf              !< index of south boundary on former subdomain
        INTEGER(iwp)       ::  nys_on_file       !< index of south boundary on former local domain
@@ -4354,5 +4422,5 @@
           CASE ( 'surfaces%var_av' )
              IF ( k == 1 )  READ ( 13 )  surfaces%var_av
-             
+
           CASE DEFAULT
 
@@ -4363,5 +4431,5 @@
 
     END SUBROUTINE surface_data_output_rrd_local
-    
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4377,5 +4445,5 @@
 
     END SUBROUTINE surface_data_output_wrd_global
-    
+
 !------------------------------------------------------------------------------!
 ! Description:
