Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3487)
+++ palm/trunk/SOURCE/Makefile	(revision 3494)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Surface output revised
+# 
+# 3474 2018-10-30 21:07:39Z kanani
 # Add virtual measurement module
 # 
@@ -1424,4 +1427,5 @@
 	salsa_mod.o \
 	surface_layer_fluxes_mod.o \
+	surface_output_mod.o \
 	time_to_string.o \
 	write_restart_data_mod.o
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 3487)
+++ palm/trunk/SOURCE/palm.f90	(revision 3494)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Last actions for surface output added
+! 
+! 3487 2018-11-05 07:18:02Z maronga
 ! Updates version to 6.0
 ! 
@@ -315,4 +318,7 @@
         ONLY:  nested_run, pmci_child_initialize, pmci_init,                   &
                pmci_modelconfiguration, pmci_parent_initialize
+               
+    USE surface_output_mod,                                                    &
+        ONLY:  surface_output_last_action
 
     USE write_restart_data_mod,                                                &
@@ -576,4 +582,8 @@
        
     ENDIF
+!
+!-- Last actions for surface output, for instantaneous and time-averaged data
+    CALL surface_output_last_action( 0 )
+    CALL surface_output_last_action( 1 )
 
 !
Index: palm/trunk/SOURCE/surface_output_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_output_mod.f90	(revision 3487)
+++ palm/trunk/SOURCE/surface_output_mod.f90	(revision 3494)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! 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
 ! 
@@ -30,12 +36,18 @@
 ! Add output variables
 ! Write data into binary file
+!
+! 3420 2018-10-24 17:30:08Z gronemeier
 ! Initial implementation from Klaus Ketelsen and Matthias Suehring
 !
-! 3420 2018-10-24 17:30:08Z gronemeier
+!
+! Authors:
+! --------
+! @author Klaus Ketelsen, Matthias Suehring
 !
 ! Description:
 ! ------------
 !> Generate output for surface data.
-!
+!> 
+!> @todo Create namelist file for post-processing tool. 
 !------------------------------------------------------------------------------!
 
@@ -46,10 +58,16 @@
    USE arrays_3d,                                                              &
        ONLY:  zu, zw
+       
    USE control_parameters,                                                     &
        ONLY:  surface_data_output
+       
    USE grid_variables,                                                         &
        ONLY: dx,dy
+       
    USE indices,                                                                &
        ONLY: nxl, nxr, nys, nyn, nzb, nzt
+       
+   USE pegrid
+       
    USE surface_mod,                                                            &
        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,                  &
@@ -58,8 +76,10 @@
    IMPLICIT NONE
 
-   TYPE surf_out                !< data structure which contains all surfaces elements of all types on subdomain
+   TYPE surf_out                      !< data structure which contains all surfaces elements of all types on subdomain
    
-      INTEGER(iwp) ::  ns       !< number of surface elements
-      INTEGER(iwp) ::  npoints  !< number of points  / vertices which define a surface element
+      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
       
       REAL(wp) ::  fillvalue = -9999.9_wp !< fillvalue for surface elements which are not defined
@@ -111,4 +131,8 @@
    END INTERFACE  surface_output_init
    
+   INTERFACE  surface_output_last_action
+      MODULE PROCEDURE surface_output_last_action
+   END INTERFACE  surface_output_last_action
+      
    INTERFACE  surface_output_parin
       MODULE PROCEDURE surface_output_parin
@@ -119,5 +143,5 @@
    PUBLIC surface_output, surface_output_averaging,                            &
           surface_output_check_parameters, surface_output_init,                &
-          surface_output_parin
+          surface_output_last_action, surface_output_parin
 !
 !--Public variables
@@ -134,17 +158,18 @@
 !------------------------------------------------------------------------------!
    SUBROUTINE surface_output_init
+   
       IMPLICIT NONE
 
-      INTEGER(iwp) ::  i     !< grid index in x-direction, also running variable for counting non-average data output
-      INTEGER(iwp) ::  ilen  !< string length
-      INTEGER(iwp) ::  j     !< grid index in y-direction, also running variable for counting average data output
-      INTEGER(iwp) ::  k     !< grid index in z-direction
-      INTEGER(iwp) ::  l     !< running index for surface-element orientation
-      INTEGER(iwp) ::  m     !< running index for surface elements
-      INTEGER(iwp) ::  n_out !< running index for number of output variables
+      INTEGER(iwp) ::  i                 !< grid index in x-direction, also running variable for counting non-average data output
+      INTEGER(iwp) ::  ilen              !< string length
+      INTEGER(iwp) ::  j                 !< grid index in y-direction, also running variable for counting average data output
+      INTEGER(iwp) ::  k                 !< grid index in z-direction
+      INTEGER(iwp) ::  l                 !< running index for surface-element orientation
+      INTEGER(iwp) ::  m                 !< running index for surface elements
+      INTEGER(iwp) ::  n_out             !< running index for number of output variables   
+      INTEGER(iwp) ::  npg               !< counter variable for all surface elements ( or polygons )
+      INTEGER(iwp) ::  point_index_count !< local counter variable for point index
       
-     
-      INTEGER(iwp) ::  npg   !< counter variable for all surface elements ( or polygons )
-      
+      INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_points_on_pe   !< array which contains the number of points on all mpi ranks
       INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) ::  point_index !< dummy array used to check where the reference points for surface polygons are located
            
@@ -187,20 +212,20 @@
 !--         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
+            IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints
+               point_index(k,j,i) = surfaces%npoints - 1
             ENDIF
-            IF ( point_index(k,j+1,i) < 0)  THEN
+            IF ( point_index(k,j,i+1) < 0 )  THEN
                surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k,j+1,i) = surfaces%npoints
+               point_index(k,j,i+1) = surfaces%npoints - 1
             ENDIF
-            IF ( point_index(k,j,i+1) < 0)  THEN
+            IF ( point_index(k,j+1,i+1) < 0 )  THEN
+               surfaces%npoints       = surfaces%npoints + 1    
+               point_index(k,j+1,i+1) = surfaces%npoints - 1
+            ENDIF
+            IF ( point_index(k,j+1,i) < 0 )  THEN
                surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k,j,i+1) = surfaces%npoints
-            ENDIF
-            IF ( point_index(k,j+1,i+1) < 0)  THEN
-               surfaces%npoints       = surfaces%npoints + 1    
-               point_index(k,j+1,i+1) = surfaces%npoints
-            ENDIF
+               point_index(k,j+1,i) = surfaces%npoints - 1
+            ENDIF           
          ENDDO
       ENDDO
@@ -210,20 +235,20 @@
          k = surf_lsm_h%k(m) + surf_lsm_h%koff
 
-         IF ( point_index(k,j,i) < 0)  THEN
+         IF ( point_index(k,j,i) < 0 )  THEN
             surfaces%npoints   = surfaces%npoints + 1    
-            point_index(k,j,i) = surfaces%npoints
+            point_index(k,j,i) = surfaces%npoints - 1
          ENDIF
-         IF ( point_index(k,j+1,i) < 0)  THEN
+         IF ( point_index(k,j,i+1) < 0 )  THEN
             surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j+1,i) = surfaces%npoints
+            point_index(k,j,i+1) = surfaces%npoints - 1
          ENDIF
-         IF ( point_index(k,j,i+1) < 0)  THEN
+         IF ( point_index(k,j+1,i+1) < 0 )  THEN
+            surfaces%npoints       = surfaces%npoints + 1    
+            point_index(k,j+1,i+1) = surfaces%npoints - 1
+         ENDIF
+         IF ( point_index(k,j+1,i) < 0 )  THEN
             surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j,i+1) = surfaces%npoints
-         ENDIF
-         IF ( point_index(k,j+1,i+1) < 0)  THEN
-            surfaces%npoints       = surfaces%npoints + 1    
-            point_index(k,j+1,i+1) = surfaces%npoints
-         ENDIF   
+            point_index(k,j+1,i) = surfaces%npoints - 1
+         ENDIF    
       ENDDO
       DO  m = 1, surf_usm_h%ns
@@ -232,19 +257,19 @@
          k = surf_usm_h%k(m) + surf_usm_h%koff
 
-         IF ( point_index(k,j,i) < 0)  THEN
+         IF ( point_index(k,j,i) < 0 )  THEN
             surfaces%npoints   = surfaces%npoints + 1    
-            point_index(k,j,i) = surfaces%npoints
+            point_index(k,j,i) = surfaces%npoints - 1
          ENDIF
-         IF ( point_index(k,j+1,i) < 0)  THEN
+         IF ( point_index(k,j,i+1) < 0 )  THEN
             surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j+1,i) = surfaces%npoints
+            point_index(k,j,i+1) = surfaces%npoints - 1
          ENDIF
-         IF ( point_index(k,j,i+1) < 0)  THEN
+         IF ( point_index(k,j+1,i+1) < 0 )  THEN
+            surfaces%npoints       = surfaces%npoints + 1    
+            point_index(k,j+1,i+1) = surfaces%npoints - 1
+         ENDIF
+         IF ( point_index(k,j+1,i) < 0 )  THEN
             surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j,i+1) = surfaces%npoints
-         ENDIF
-         IF ( point_index(k,j+1,i+1) < 0)  THEN
-            surfaces%npoints       = surfaces%npoints + 1    
-            point_index(k,j+1,i+1) = surfaces%npoints
+            point_index(k,j+1,i) = surfaces%npoints - 1
          ENDIF     
       ENDDO
@@ -258,124 +283,125 @@
 !--         identical to the reference j-index outside the grid box. 
 !--         Equivalent for east-facing surfaces and i-index.
-            i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 2 )
-            j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 0 )
+            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
+            IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints
+               point_index(k,j,i) = surfaces%npoints - 1
             ENDIF 
 !
+!--         Upper / lower right index for north- and south-facing surfaces
+            IF ( l == 0  .OR.  l == 1 )  THEN
+               IF ( point_index(k,j,i+1) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j,i+1) = surfaces%npoints - 1
+               ENDIF 
+               IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j,i+1) = surfaces%npoints - 1
+               ENDIF
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+               IF ( point_index(k,j+1,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j+1,i) = surfaces%npoints - 1
+               ENDIF 
+               IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j+1,i) = surfaces%npoints - 1
+               ENDIF
+            ENDIF
+!
 !--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0)  THEN
+            IF ( point_index(k+1,j,i) < 0 )  THEN
                surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints
+               point_index(k+1,j,i) = surfaces%npoints - 1
             ENDIF
+         ENDDO
+         DO  m = 1, surf_lsm_v(l)%ns
+            i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
+            j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
+            k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
+!
+!--         Lower left /front vertex
+            IF ( point_index(k,j,i) < 0 )  THEN
+               surfaces%npoints   = surfaces%npoints + 1    
+               point_index(k,j,i) = surfaces%npoints - 1
+            ENDIF 
 !
 !--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .AND.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0)  THEN
+            IF ( l == 0  .OR.  l == 1 )  THEN
+               IF ( point_index(k,j,i+1) < 0 )  THEN
                   surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints
+                  point_index(k,j,i+1) = surfaces%npoints - 1
                ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0)  THEN
+               IF ( point_index(k+1,j,i+1) < 0 )  THEN
                   surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints
+                  point_index(k+1,j,i+1) = surfaces%npoints - 1
                ENDIF
 !
 !--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .AND.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0)  THEN
+            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+               IF ( point_index(k,j+1,i) < 0 )  THEN
                   surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints
+                  point_index(k,j+1,i) = surfaces%npoints - 1
                ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0)  THEN
+               IF ( point_index(k+1,j+1,i) < 0 )  THEN
                   surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints
-               ENDIF
-            ENDIF            
-            
-         ENDDO
-         DO  m = 1, surf_lsm_v(l)%ns
-            i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 2 )
-            j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 0 )
-            k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
+                  point_index(k+1,j+1,i) = surfaces%npoints - 1
+               ENDIF
+            ENDIF
+!
+!--         Upper left / front vertex
+            IF ( point_index(k+1,j,i) < 0 )  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k+1,j,i) = surfaces%npoints - 1
+            ENDIF
+         ENDDO
+
+         DO  m = 1, surf_usm_v(l)%ns
+            i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
+            j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
+            k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
 !
 !--         Lower left /front vertex
-            IF ( point_index(k,j,i) < 0)  THEN
+            IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints
+               point_index(k,j,i) = surfaces%npoints - 1
             ENDIF 
 !
+!--         Upper / lower right index for north- and south-facing surfaces
+            IF ( l == 0  .OR.  l == 1 )  THEN
+               IF ( point_index(k,j,i+1) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j,i+1) = surfaces%npoints - 1
+               ENDIF 
+               IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j,i+1) = surfaces%npoints - 1
+               ENDIF
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+               IF ( point_index(k,j+1,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j+1,i) = surfaces%npoints - 1
+               ENDIF 
+               IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j+1,i) = surfaces%npoints - 1
+               ENDIF
+            ENDIF
+!
 !--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0)  THEN
+            IF ( point_index(k+1,j,i) < 0 )  THEN
                surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints
+               point_index(k+1,j,i) = surfaces%npoints - 1
             ENDIF
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .AND.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .AND.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints
-               ENDIF
-            ENDIF       
-         ENDDO
-         DO  m = 1, surf_usm_v(l)%ns
-            i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 2 )
-            j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 0 )
-            k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
-!
-!--         Lower left /front vertex
-            IF ( point_index(k,j,i) < 0)  THEN
-               surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints
-            ENDIF 
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0)  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints
-            ENDIF
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .AND.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .AND.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints
-               ENDIF
-            ENDIF    
-         ENDDO
+         ENDDO
+
       ENDDO
 !
@@ -384,12 +410,34 @@
 !--   of points (vertices) which define the polygons can be larger.
       ALLOCATE( surfaces%points(3,1:surfaces%npoints) )
-      ALLOCATE( surfaces%polygons(5,1:surfaces%ns)    )      
+      ALLOCATE( surfaces%polygons(5,1:surfaces%ns)    )
+!
+!--   Note, PARAVIEW expects consecutively ordered points which can be
+!--   unambiguously identified. 
+!--   Hence, all PEs know where they should start counting, depending on the
+!--   number of points on the other PE's with lower MPI rank. 
+#if defined( __parallel )
+      CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER,                    &
+                          num_points_on_pe, 1, MPI_INTEGER, comm2d, ierr  )
+#else
+      num_points_on_pe = surfaces%npoints
+#endif
+
 !
 !--   After the number of vertices is counted, repeat the loops and define the
-!--   vertices. 
-!--   Start with the horizontal default surfaces
+!--   vertices. Start with the horizontal default surfaces
+!--   First, however, determine the offset where couting of points should be
+!--   started, which is the sum of points of all PE's with lower MPI rank.
+      i                 = 0
+      point_index_count = 0
+      DO WHILE ( i < myid  .AND.  i <= SIZE( num_points_on_pe ) )
+         point_index_count = point_index_count + num_points_on_pe(i)
+         i                 = i + 1
+      ENDDO
+          
+      
       surfaces%npoints = 0
       point_index      = -1
       npg              = 0
+      
       DO  l = 0, 1
          DO  m = 1, surf_def_h(0)%ns
@@ -402,39 +450,43 @@
 !--         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
+            IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints
+               point_index(k,j,i) = point_index_count
+               point_index_count  = point_index_count + 1               
                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                surfaces%points(3,surfaces%npoints) = zw(k)
             ENDIF
-            IF ( point_index(k,j+1,i) < 0)  THEN
+            IF ( point_index(k,j,i+1) < 0 )  THEN
                surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k,j+1,i) = surfaces%npoints
+               point_index(k,j,i+1) = point_index_count
+               point_index_count    = point_index_count + 1 
+               surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+               surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+               surfaces%points(3,surfaces%npoints) = zw(k)
+            ENDIF
+            IF ( point_index(k,j+1,i+1) < 0 )  THEN
+               surfaces%npoints       = surfaces%npoints + 1    
+               point_index(k,j+1,i+1) = point_index_count
+               point_index_count      = point_index_count + 1 
+               surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+               surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+               surfaces%points(3,surfaces%npoints) = zw(k)
+            ENDIF
+            IF ( point_index(k,j+1,i) < 0 )  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k,j+1,i) = point_index_count
+               point_index_count    = point_index_count + 1 
                surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
                surfaces%points(3,surfaces%npoints) = zw(k)
             ENDIF
-            IF ( point_index(k,j,i+1) < 0)  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k,j,i+1) = surfaces%npoints
-               surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k)
-            ENDIF
-            IF ( point_index(k,j+1,i+1) < 0)  THEN
-               surfaces%npoints       = surfaces%npoints + 1    
-               point_index(k,j+1,i+1) = surfaces%npoints
-               surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k)
-            ENDIF
             
-            npg = npg + 1
+            npg                        = npg + 1
             surfaces%polygons(1,npg)   = 4
             surfaces%polygons(2,npg)   = point_index(k,j,i)
-            surfaces%polygons(3,npg)   = point_index(k,j+1,i)
+            surfaces%polygons(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,i+1)  
+            surfaces%polygons(5,npg)   = point_index(k,j+1,i)
          ENDDO
       ENDDO
@@ -443,82 +495,91 @@
          j = surf_lsm_h%j(m) + surf_lsm_h%joff
          k = surf_lsm_h%k(m) + surf_lsm_h%koff
-
-         IF ( point_index(k,j,i) < 0)  THEN
+         IF ( point_index(k,j,i) < 0 )  THEN
             surfaces%npoints   = surfaces%npoints + 1    
-            point_index(k,j,i) = surfaces%npoints
+            point_index(k,j,i) = point_index_count
+            point_index_count  = point_index_count + 1 
             surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
             surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
             surfaces%points(3,surfaces%npoints) = zw(k)
          ENDIF
-         IF ( point_index(k,j+1,i) < 0)  THEN
+         IF ( point_index(k,j,i+1) < 0 )  THEN
             surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j+1,i) = surfaces%npoints
+            point_index(k,j,i+1) = point_index_count
+            point_index_count    = point_index_count + 1 
+            surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+            surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+            surfaces%points(3,surfaces%npoints) = zw(k)
+         ENDIF
+         IF ( point_index(k,j+1,i+1) < 0 )  THEN
+            surfaces%npoints       = surfaces%npoints + 1    
+            point_index(k,j+1,i+1) = point_index_count
+            point_index_count      = point_index_count + 1 
+            surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+            surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+            surfaces%points(3,surfaces%npoints) = zw(k)
+         ENDIF
+         IF ( point_index(k,j+1,i) < 0 )  THEN
+            surfaces%npoints     = surfaces%npoints + 1    
+            point_index(k,j+1,i) = point_index_count
+            point_index_count    = point_index_count + 1 
             surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
             surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
             surfaces%points(3,surfaces%npoints) = zw(k)
          ENDIF
-         IF ( point_index(k,j,i+1) < 0)  THEN
+         
+         npg                        = npg + 1
+         surfaces%polygons(1,npg)   = 4
+         surfaces%polygons(2,npg)   = point_index(k,j,i)
+         surfaces%polygons(3,npg)   = point_index(k,j,i+1)
+         surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
+         surfaces%polygons(5,npg)   = point_index(k,j+1,i) 
+      ENDDO
+
+      DO  m = 1, surf_usm_h%ns
+         i = surf_usm_h%i(m) + surf_usm_h%ioff
+         j = surf_usm_h%j(m) + surf_usm_h%joff
+         k = surf_usm_h%k(m) + surf_usm_h%koff
+
+         IF ( point_index(k,j,i) < 0 )  THEN
+            surfaces%npoints   = surfaces%npoints + 1    
+            point_index(k,j,i) = point_index_count
+            point_index_count  = point_index_count + 1 
+            surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+            surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
+            surfaces%points(3,surfaces%npoints) = zw(k)
+         ENDIF
+         IF ( point_index(k,j,i+1) < 0 )  THEN
             surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j,i+1) = surfaces%npoints
+            point_index(k,j,i+1) = point_index_count
+            point_index_count    = point_index_count + 1 
             surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
             surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
             surfaces%points(3,surfaces%npoints) = zw(k)
          ENDIF
-         IF ( point_index(k,j+1,i+1) < 0)  THEN
+         IF ( point_index(k,j+1,i+1) < 0 )  THEN
             surfaces%npoints       = surfaces%npoints + 1    
-            point_index(k,j+1,i+1) = surfaces%npoints
+            point_index(k,j+1,i+1) = point_index_count
+            point_index_count      = point_index_count + 1 
             surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
             surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
             surfaces%points(3,surfaces%npoints) = zw(k)
          ENDIF
-         
-         npg = npg + 1
-         surfaces%polygons(1,npg)   = 4
-         surfaces%polygons(2,npg)   = point_index(k,j,i)
-         surfaces%polygons(3,npg)   = point_index(k,j+1,i)
-         surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
-         surfaces%polygons(5,npg)   = point_index(k,j,i+1)  
-      ENDDO
-      DO  m = 1, surf_usm_h%ns
-         i = surf_usm_h%i(m) + surf_usm_h%ioff
-         j = surf_usm_h%j(m) + surf_usm_h%joff
-         k = surf_usm_h%k(m) + surf_usm_h%koff
-
-         IF ( point_index(k,j,i) < 0)  THEN
-            surfaces%npoints   = surfaces%npoints + 1    
-            point_index(k,j,i) = surfaces%npoints
-            surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         IF ( point_index(k,j+1,i) < 0)  THEN
+         IF ( point_index(k,j+1,i) < 0 )  THEN
             surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j+1,i) = surfaces%npoints
+            point_index(k,j+1,i) = point_index_count
+            point_index_count    = point_index_count + 1 
             surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
             surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
             surfaces%points(3,surfaces%npoints) = zw(k)
          ENDIF
-         IF ( point_index(k,j,i+1) < 0)  THEN
-            surfaces%npoints     = surfaces%npoints + 1    
-            point_index(k,j,i+1) = surfaces%npoints
-            surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
-         IF ( point_index(k,j+1,i+1) < 0)  THEN
-            surfaces%npoints       = surfaces%npoints + 1    
-            point_index(k,j+1,i+1) = surfaces%npoints
-            surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-            surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-            surfaces%points(3,surfaces%npoints) = zw(k)
-         ENDIF
          
-         npg = npg + 1
+         npg                        = npg + 1
          surfaces%polygons(1,npg)   = 4
          surfaces%polygons(2,npg)   = point_index(k,j,i)
-         surfaces%polygons(3,npg)   = point_index(k,j+1,i)
+         surfaces%polygons(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,i+1)    
+         surfaces%polygons(5,npg)   = point_index(k,j+1,i)    
       ENDDO
+
       DO  l = 0, 3
          DO  m = 1, surf_def_v(l)%ns
@@ -528,63 +589,69 @@
 !--         identical to the reference j-index outside the grid box. 
 !--         Equivalent for east-facing surfaces and i-index.
-            i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 2 )
-            j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 0 )
+            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
+            IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints
+               point_index(k,j,i) = point_index_count
+               point_index_count  = point_index_count + 1 
+               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
+               surfaces%points(3,surfaces%npoints) = zw(k-1)
+            ENDIF 
+!
+!--         Upper / lower right index for north- and south-facing surfaces
+            IF ( l == 0  .OR.  l == 1 )  THEN
+               IF ( point_index(k,j,i+1) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j,i+1) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k-1)
+               ENDIF 
+               IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j,i+1) = point_index_count
+                  point_index_count      = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+               IF ( point_index(k,j+1,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j+1,i) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k-1)
+               ENDIF 
+               IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j+1,i) = point_index_count
+                  point_index_count      = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
+            ENDIF
+!
+!--         Upper left / front vertex
+            IF ( point_index(k+1,j,i) < 0 )  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k+1,j,i) = point_index_count
+               point_index_count    = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                surfaces%points(3,surfaces%npoints) = zw(k)
-            ENDIF 
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0)  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints
-               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k+1)
             ENDIF
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .AND.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k+1)
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .AND.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k+1)
-               ENDIF
-            ENDIF          
             
             npg = npg + 1            
-            IF ( l == 0  .AND.  l == 1 )  THEN
+            IF ( l == 0  .OR.  l == 1 )  THEN
                surfaces%polygons(1,npg)   = 4
                surfaces%polygons(2,npg)   = point_index(k,j,i)
@@ -601,64 +668,71 @@
             
          ENDDO
+
          DO  m = 1, surf_lsm_v(l)%ns
-            i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 2 )
-            j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 0 )
+            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
+            IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints
+               point_index(k,j,i) = point_index_count
+               point_index_count  = point_index_count + 1
+               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
+               surfaces%points(3,surfaces%npoints) = zw(k-1)
+            ENDIF 
+!
+!--         Upper / lower right index for north- and south-facing surfaces
+            IF ( l == 0  .OR.  l == 1 )  THEN
+               IF ( point_index(k,j,i+1) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j,i+1) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k-1)
+               ENDIF 
+               IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j,i+1) = point_index_count
+                  point_index_count      = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+               IF ( point_index(k,j+1,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j+1,i) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k-1)
+               ENDIF 
+               IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j+1,i) = point_index_count
+                  point_index_count      = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
+            ENDIF
+!
+!--         Upper left / front vertex
+            IF ( point_index(k+1,j,i) < 0 )  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k+1,j,i) = point_index_count
+               point_index_count    = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                surfaces%points(3,surfaces%npoints) = zw(k)
-            ENDIF 
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0)  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints
-               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k+1)
-            ENDIF
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .AND.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k+1)
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .AND.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k+1)
-               ENDIF
             ENDIF
             
             npg = npg + 1            
-            IF ( l == 0  .AND.  l == 1 )  THEN
+            IF ( l == 0  .OR.  l == 1 )  THEN
                surfaces%polygons(1,npg)   = 4
                surfaces%polygons(2,npg)   = point_index(k,j,i)
@@ -675,63 +749,69 @@
          ENDDO
          DO  m = 1, surf_usm_v(l)%ns
-            i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 2 )
-            j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 0 )
+            i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
+            j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
             k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
 !
 !--         Lower left /front vertex
-            IF ( point_index(k,j,i) < 0)  THEN
+            IF ( point_index(k,j,i) < 0 )  THEN
                surfaces%npoints   = surfaces%npoints + 1    
-               point_index(k,j,i) = surfaces%npoints
+               point_index(k,j,i) = point_index_count
+               point_index_count  = point_index_count + 1
+               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
+               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
+               surfaces%points(3,surfaces%npoints) = zw(k-1)
+            ENDIF 
+!
+!--         Upper / lower right index for north- and south-facing surfaces
+            IF ( l == 0  .OR.  l == 1 )  THEN
+               IF ( point_index(k,j,i+1) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j,i+1) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k-1)
+               ENDIF 
+               IF ( point_index(k+1,j,i+1) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j,i+1) = point_index_count
+                  point_index_count      = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
+!
+!--         Upper / lower front index for east- and west-facing surfaces
+            ELSEIF ( l == 2  .OR.  l == 3 )  THEN
+               IF ( point_index(k,j+1,i) < 0 )  THEN
+                  surfaces%npoints     = surfaces%npoints + 1    
+                  point_index(k,j+1,i) = point_index_count
+                  point_index_count    = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k-1)
+               ENDIF 
+               IF ( point_index(k+1,j+1,i) < 0 )  THEN
+                  surfaces%npoints       = surfaces%npoints + 1    
+                  point_index(k+1,j+1,i) = point_index_count
+                  point_index_count      = point_index_count + 1
+                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
+                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
+                  surfaces%points(3,surfaces%npoints) = zw(k)
+               ENDIF
+            ENDIF
+!
+!--         Upper left / front vertex
+            IF ( point_index(k+1,j,i) < 0 )  THEN
+               surfaces%npoints     = surfaces%npoints + 1    
+               point_index(k+1,j,i) = point_index_count
+               point_index_count    = point_index_count + 1
                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
                surfaces%points(3,surfaces%npoints) = zw(k)
-            ENDIF 
-!
-!--         Upper left / front vertex
-            IF ( point_index(k+1,j,i) < 0)  THEN
-               surfaces%npoints     = surfaces%npoints + 1    
-               point_index(k+1,j,i) = surfaces%npoints
-               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
-               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
-               surfaces%points(3,surfaces%npoints) = zw(k+1)
-            ENDIF
-!
-!--         Upper / lower right index for north- and south-facing surfaces
-            IF ( l == 0  .AND.  l == 1 )  THEN
-               IF ( point_index(k,j,i+1) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j,i+1) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF 
-               IF ( point_index(k+1,j,i+1) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j,i+1) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k+1)
-               ENDIF
-!
-!--         Upper / lower front index for east- and west-facing surfaces
-            ELSEIF ( l == 2  .AND.  l == 3 )  THEN
-               IF ( point_index(k,j+1,i) < 0)  THEN
-                  surfaces%npoints     = surfaces%npoints + 1    
-                  point_index(k,j+1,i) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k)
-               ENDIF 
-               IF ( point_index(k+1,j+1,i) < 0)  THEN
-                  surfaces%npoints       = surfaces%npoints + 1    
-                  point_index(k+1,j+1,i) = surfaces%npoints
-                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
-                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
-                  surfaces%points(3,surfaces%npoints) = zw(k+1)
-               ENDIF
             ENDIF
             
             npg = npg + 1            
-            IF ( l == 0  .AND.  l == 1 )  THEN
+            IF ( l == 0  .OR.  l == 1 )  THEN
                surfaces%polygons(1,npg)   = 4
                surfaces%polygons(2,npg)   = point_index(k,j,i)
@@ -747,4 +827,5 @@
             ENDIF
          ENDDO
+
       ENDDO
 !
@@ -752,6 +833,17 @@
       DEALLOCATE ( point_index )
 !
-!--   Now, vertices as well as polygon data can be written to NetCDF files
-!--   ... 
+!--   Sum-up total number of surface elements and vertices on domain. This
+!--   will be needed for post-processing. 
+      surfaces%npoints_total = 0
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1,        &
+                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
+       CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1,                  &
+                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
+#else
+       surfaces%npoints_total = surfaces%npoints
+       surfaces%ns_total      = surfaces%ns
+#endif
+
    END SUBROUTINE surface_output_init
    
@@ -795,5 +887,7 @@
             IF ( i == io_group )  THEN
                WRITE ( 25+av )  surfaces%npoints
+               WRITE ( 25+av )  surfaces%npoints_total
                WRITE ( 25+av )  surfaces%ns
+               WRITE ( 25+av )  surfaces%ns_total
                WRITE ( 25+av )  surfaces%points
                WRITE ( 25+av )  surfaces%polygons
@@ -805,17 +899,4 @@
          ENDDO
       ENDIF
-!
-!--   Write time coordinate
-      DO  i = 0, io_blocks-1
-         IF ( i == io_group )  THEN
-            trimvar = 'time'
-            WRITE ( 25+av )  trimvar
-            WRITE ( 25+av )  simulated_time
-         ENDIF
-#if defined( __parallel )
-         CALL MPI_BARRIER( comm2d, ierr )
-#endif
-      ENDDO
-
 
       n_out = 0
@@ -835,19 +916,19 @@
                IF ( av == 0 )  THEN
                   CALL surface_output_collect( surf_def_h(0)%us,               &
-                                             surf_def_h(1)%us,                 &
-                                             surf_lsm_h%us,                    &
-                                             surf_usm_h%us,                    &
-                                             surf_def_v(0)%us,                 &
-                                             surf_lsm_v(0)%us,                 &
-                                             surf_usm_v(0)%us,                 &
-                                             surf_def_v(1)%us,                 &
-                                             surf_lsm_v(1)%us,                 &
-                                             surf_usm_v(1)%us,                 &
-                                             surf_def_v(2)%us,                 &
-                                             surf_lsm_v(2)%us,                 &
-                                             surf_usm_v(2)%us,                 &
-                                             surf_def_v(3)%us,                 &
-                                             surf_lsm_v(3)%us,                 &
-                                             surf_usm_v(3)%us ) 
+                                               surf_def_h(1)%us,               &
+                                               surf_lsm_h%us,                  &
+                                               surf_usm_h%us,                  &
+                                               surf_def_v(0)%us,               &
+                                               surf_lsm_v(0)%us,               &
+                                               surf_usm_v(0)%us,               &
+                                               surf_def_v(1)%us,               &
+                                               surf_lsm_v(1)%us,               &
+                                               surf_usm_v(1)%us,               &
+                                               surf_def_v(2)%us,               &
+                                               surf_lsm_v(2)%us,               &
+                                               surf_usm_v(2)%us,               &
+                                               surf_def_v(3)%us,               &
+                                               surf_lsm_v(3)%us,               &
+                                               surf_usm_v(3)%us ) 
                ELSE
 !
@@ -1945,5 +2026,9 @@
          DO  i = 0, io_blocks-1
             IF ( i == io_group )  THEN
-               WRITE ( 25+av )  trimvar
+               WRITE ( 25+av )  LEN_TRIM( 'time' )
+               WRITE ( 25+av )  'time'
+               WRITE ( 25+av )  simulated_time
+               WRITE ( 25+av )  LEN_TRIM( trimvar )
+               WRITE ( 25+av )  TRIM( trimvar )
                WRITE ( 25+av )  surfaces%var_out
             ENDIF
@@ -1956,6 +2041,6 @@
       
 !
-!--      If averaged output was written to NetCDF file, set the counter to zero
-         IF ( av == 1 )  average_count_surf = 0
+!--   If averaged output was written to NetCDF file, set the counter to zero
+      IF ( av == 1 )  average_count_surf = 0
              
    END SUBROUTINE surface_output
@@ -2679,12 +2764,10 @@
 !> Sum-up the surface data for average output variables.  
 !------------------------------------------------------------------------------!   
-   SUBROUTINE surface_output_sum_up(  var_def_h0, var_def_h1,                  &
-                                      var_lsm_h, var_usm_h,                    &
-                                      var_def_v0, var_def_v1, var_def_v2,      &
-                                      var_def_v3,                              &
-                                      var_lsm_v0, var_lsm_v1, var_lsm_v2,      &
-                                      var_lsm_v3,                              &
-                                      var_usm_v0, var_usm_v1, var_usm_v2,      &
-                                      var_usm_v3, n_out )
+   SUBROUTINE surface_output_sum_up( var_def_h0, var_def_h1,                   &
+                                     var_lsm_h,  var_usm_h,                    &
+                                     var_def_v0, var_lsm_v0, var_usm_v0,       &
+                                     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, n_out )
    
       IMPLICIT NONE
@@ -2886,11 +2969,9 @@
 !------------------------------------------------------------------------------!
    SUBROUTINE surface_output_collect( var_def_h0, var_def_h1,                  &
-                                      var_lsm_h, var_usm_h,                    &
-                                      var_def_v0, var_def_v1, var_def_v2,      &
-                                      var_def_v3,                              &
-                                      var_lsm_v0, var_lsm_v1, var_lsm_v2,      &
-                                      var_lsm_v3,                              &
-                                      var_usm_v0, var_usm_v1, var_usm_v2,      &
-                                      var_usm_v3 )
+                                      var_lsm_h,  var_usm_h,                   &
+                                      var_def_v0, var_lsm_v0, var_usm_v0,      &
+                                      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 )                                               
    
       IMPLICIT NONE
@@ -3246,4 +3327,43 @@
 
     END SUBROUTINE surface_output_check_parameters
+    
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Last action.
+!------------------------------------------------------------------------------!
+    SUBROUTINE surface_output_last_action( av )
+  
+      USE control_parameters,                                                  &
+          ONLY:  io_blocks, io_group
+
+      USE pegrid,                                                              &
+          ONLY:  comm2d, ierr
+
+      IMPLICIT NONE
+      
+      INTEGER(iwp) ::  av     !< id indicating average or non-average data output
+      INTEGER(iwp) ::  i      !< loop index
+
+!
+!--   Return, if nothing to output
+      IF ( dosurf_no(av) == 0 )  RETURN
+!
+!--   Open files
+      CALL check_open( 25+av )
+!
+!--   Write time coordinate
+      DO  i = 0, io_blocks-1
+         IF ( i == io_group )  THEN
+            WRITE ( 25+av )  LEN_TRIM( 'END' )
+            WRITE ( 25+av )  'END'
+         ENDIF
+#if defined( __parallel )
+         CALL MPI_BARRIER( comm2d, ierr )
+#endif
+      ENDDO    
+
+    END SUBROUTINE surface_output_last_action    
 
 !--Private Subroutines (Only called from tinside his module)
Index: palm/trunk/SOURCE/virtual_measurement_mod.f90
===================================================================
--- palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 3487)
+++ palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 3494)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfixing
+! 
+! 3473 2018-10-30 20:50:15Z suehring
 ! Initial revision
 ! 
@@ -31,5 +34,5 @@
 ! Authors:
 ! --------
-! @author Matthias Suehring
+! @author Matthias Suehring and Klaus Ketelsen
 ! 
 !
@@ -554,14 +557,9 @@
 !--          lower-left origin of the model domain, which is the difference 
 !--          betwen UTM coordinates.
-!              write(9,*) l, "before Eutm", e_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_x
-!              write(9,*) l, "before Nutm", n_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_y
-             
+            
              e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t))         &
                                          - init_model%origin_x
              n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t))         &
                                          - init_model%origin_y
-!              write(9,*) l, "Eutm", e_utm(t,1:vmea(l)%dim_t(t)) 
-!              write(9,*) l, "Nutm", n_utm(t,1:vmea(l)%dim_t(t)) 
-!              write(9,*)
 !
 !--          Compute grid indices relative to origin and check if these are
@@ -585,5 +583,4 @@
                 ks = k_prev
                 IF ( on_pe )  THEN
-!                    write(9,*) " I am onpe" 
                    ksurf = get_topography_top_index_ji( js, is, 's' )
                    ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
@@ -655,5 +652,4 @@
                    vmea(l)%j(ns)  = js
                    vmea(l)%k(ns)  = ks
-!                    write(9,*) TRIM( vmea(l)%feature_type), l, "ns", ns, "ijk", vmea(l)%i(ns), vmea(l)%j(ns), vmea(l)%k(ns)
                 ENDIF
 !
@@ -713,6 +709,4 @@
      INTEGER(iwp) ::  mm, j, i
 
-!      write(9,*) "sampling"
-!      flush(9)
 !
 !--  Loop over all stations. For each possible variable loop over all 
@@ -723,7 +717,7 @@
 !--     the same indices for scalar and velocity components are used.
 !--     ToDo: Revise this later.
-!         DO  m = 1, vmea(l)%ns
-!            j = vmea(l)%j(m)
-!            i = vmea(l)%i(m)
+        DO  m = 1, vmea(l)%ns
+           j = vmea(l)%j(m)
+           i = vmea(l)%i(m)
 !            
 !            IF ( i >= nxl  .AND.  i <= nxr  .AND.                               &
@@ -731,7 +725,4 @@
 !               IF ( surf_def_h(0)%start_index(j,i) <= &
 !                    surf_def_h(0)%end_index(j,i) )  THEN
-!                 
-!                  write(9,*) "sampled"
-!                  flush(9)
 !                  mm = surf_def_h(0)%end_index(j,i)
 !               
@@ -763,5 +754,5 @@
 !               ENDDO
 !            ENDIF
-!         ENDDO
+        ENDDO
 
      ENDDO
