Index: /palm/trunk/SOURCE/surface_data_output_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_data_output_mod.f90	(revision 4599)
+++ /palm/trunk/SOURCE/surface_data_output_mod.f90	(revision 4600)
@@ -20,9 +20,15 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! - Change: adjustmens for mpi-io - surface data is transformed to a 2D-based surface array
+!   before writing.
+! - Bugfix in counting of surface elements
+! - Bugfix in data-output of averaged surface data in case of restarts
+! 
+! 4577 2020-06-25 09:53:58Z raasch
 ! File re-formatted to follow the PALM coding standard
 !
@@ -178,7 +184,10 @@
 
    USE restart_data_mpi_io_mod,                                                                    &
-       ONLY:  rd_mpi_io_check_array,                                                               &
-              rrd_mpi_io,                                                                          &
-              wrd_mpi_io
+       ONLY:  rrd_mpi_io,                                                                          &
+              rd_mpi_io_check_array,                                                               &
+              rrd_mpi_io_surface,                                                                  &
+              rd_mpi_io_surface_filetypes,                                                         &
+              wrd_mpi_io,                                                                          &
+              wrd_mpi_io_surface
 
    USE surface_mod,                                                                                &
@@ -213,6 +222,6 @@
       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
+      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_out   !< output variable
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_av    !< variable used for averaging
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  points    !< points  / vertices of a surface element
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  polygons  !< polygon data of a surface element
@@ -398,7 +407,7 @@
 
 #if defined( __netcdf4_parallel )
-    CHARACTER (LEN=100)  ::  filename            !< name of output file
-    CHARACTER (LEN=80)   ::   time_average_text  !< string written to file attribute time_avg
-    CHARACTER (LEN=4000) ::   var_list           !< list of variables written to NetCDF file
+    CHARACTER (LEN=100)  ::  filename           !< name of output file
+    CHARACTER (LEN=80)   ::  time_average_text  !< string written to file attribute time_avg
+    CHARACTER (LEN=4000) ::  var_list           !< list of variables written to NetCDF file
 
     INTEGER(iwp) ::  av                 !< flag for averaged (=1) and non-averaged (=0) data
@@ -436,10 +445,10 @@
        surfaces%npoints = 0
        DO  l = 0, 1
-          DO  m = 1, surf_def_h(0)%ns
+          DO  m = 1, surf_def_h(l)%ns
 !
 !--          Determine the indices of the respective grid cell inside the topography
-             i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
-             j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
-             k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
+             i = surf_def_h(l)%i(m) + surf_def_h(l)%ioff
+             j = surf_def_h(l)%j(m) + surf_def_h(l)%joff
+             k = surf_def_h(l)%k(m) + surf_def_h(l)%koff
 !
 !--          Check if the vertices that define the surface element are already defined, if not,
@@ -4579,4 +4588,6 @@
        CASE ( 'average_count_surf' )
           READ ( 13 )  average_count_surf
+       CASE ( 'time_dosurf_av' )
+          READ ( 13 )  time_dosurf_av
 
        CASE DEFAULT
@@ -4598,4 +4609,5 @@
 
     CALL rrd_mpi_io( 'average_count_surf', average_count_surf )
+    CALL rrd_mpi_io( 'time_dosurf_av', time_dosurf_av )
 
  END SUBROUTINE surface_data_output_rrd_global_mpi
@@ -4645,14 +4657,111 @@
     IMPLICIT NONE
 
-    LOGICAL ::  array_found  !<
-
-
-    CALL rd_mpi_io_check_array( 'surfaces%var_av' , found = array_found )
-
-!> does not work this way: surface%var_av has non-standard dimensions
-!    IF ( array_found )  THEN
-!       IF ( .NOT. ALLOCATED( surfaces%var_av ) )  ALLOCATE( ....... )
-!       CALL rrd_mpi_io( 'surfaces%var_av', surfaces%var_av )
-!    ENDIF
+    CHARACTER(LEN=3) ::  dum !< dummy string to create output-variable name
+
+    INTEGER(iwp) ::  i  !< grid index in x-direction
+    INTEGER(iwp) ::  j  !< grid index in y-direction
+    INTEGER(iwp) ::  l  !< running index surface orientation
+    INTEGER(iwp) ::  m  !< running index surface elements
+    INTEGER(iwp) ::  n  !< counting variable
+    INTEGER(iwp) ::  nv !< running index over number of variables
+
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  end_index           !< end index of surface data at (j,i)
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  global_start_index  !< index array for surface data (MPI-IO)
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  num_surf            !< number of surface data at (j,i)
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  start_index         !< start index of surface data at (j,i)
+
+    LOGICAL ::  array_found  !< flag indicating whether variable is in restart data or not
+    LOGICAL ::  ldum         !< dummy variable
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surf_in !< input array in expected restart format
+
+!
+!-- Note, surface data which is written to file is organized in a different way than
+!-- the output surface data. The output surface data is a concatenated array of the 
+!-- different surface types and orientations, while the mpi-io expects surface data that
+!-- is consecutive in terms of start- and end-index, i.e. organized along the (j,i)
+!-- grid index. Hence, data need to be tranformed back to the output surface data.
+    ALLOCATE( end_index(nys:nyn,nxl:nxr)          )
+    ALLOCATE( num_surf(nys:nyn,nxl:nxr)           )
+    ALLOCATE( start_index(nys:nyn,nxl:nxr)        )
+    ALLOCATE( global_start_index(nys:nyn,nxl:nxr) )
+
+    ALLOCATE( surf_in(1:surfaces%ns) )
+
+    CALL rd_mpi_io_check_array( 'surfaces%start_index', found = array_found )
+    IF ( array_found )  CALL rrd_mpi_io( 'surfaces%start_index', start_index )
+
+    CALL rd_mpi_io_check_array( 'surfaces%end_index', found = array_found )
+    IF ( array_found )  CALL rrd_mpi_io( 'surfaces%end_index', end_index )
+
+    CALL rd_mpi_io_check_array( 'surfaces%global_start_index', found = array_found )
+    IF ( array_found )  CALL rrd_mpi_io( 'surfaces%global_start_index', global_start_index )
+
+    CALL rd_mpi_io_surface_filetypes( start_index, end_index, ldum, global_start_index )
+
+    DO  nv = 1, dosurf_no(1)
+       IF ( nv < 10                     )  WRITE( dum, '(I1)') nv
+       IF ( nv < 100   .AND.  nv >= 10  )  WRITE( dum, '(I2)') nv
+       IF ( nv < 1000  .AND.  nv >= 100 )  WRITE( dum, '(I3)') nv
+
+       CALL rd_mpi_io_check_array( 'surfaces%var_av' // TRIM( dum ), found = array_found )
+
+       IF ( array_found )  THEN
+
+          CALL rrd_mpi_io_surface( 'surfaces%var_av' // TRIM(dum), surf_in )
+!
+!--       Write temporary input variable back to surface-output data array.
+          n = 0
+          num_surf = 0
+          DO  l = 0, 1
+             DO  m = 1, surf_def_h(l)%ns
+                i = surf_def_h(l)%i(m)
+                j = surf_def_h(l)%j(m)
+                n                     = n + 1
+                surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
+                num_surf(j,i)         = num_surf(j,i) + 1
+             ENDDO
+          ENDDO
+          DO  m = 1, surf_lsm_h%ns
+             i = surf_lsm_h%i(m)
+             j = surf_lsm_h%j(m)
+             n                     = n + 1
+             surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
+             num_surf(j,i)         = num_surf(j,i) + 1
+          ENDDO
+          DO  m = 1, surf_usm_h%ns 
+             i = surf_usm_h%i(m)
+             j = surf_usm_h%j(m)
+             n                     = n + 1
+             surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
+             num_surf(j,i)         = num_surf(j,i) + 1
+          ENDDO
+
+          DO  l = 0, 3
+             DO  m = 1, surf_def_v(l)%ns
+                i = surf_def_v(l)%i(m)
+                j = surf_def_v(l)%j(m)
+                n                     = n + 1
+                surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
+                num_surf(j,i)         = num_surf(j,i) + 1
+             ENDDO
+             DO  m = 1, surf_lsm_v(l)%ns
+                i = surf_lsm_v(l)%i(m)
+                j = surf_lsm_v(l)%j(m)
+                n                     = n + 1
+                surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
+                num_surf(j,i)         = num_surf(j,i) + 1
+             ENDDO
+             DO  m = 1, surf_usm_v(l)%ns
+                i = surf_usm_v(l)%i(m)
+                j = surf_usm_v(l)%j(m)
+                n                     = n + 1
+                surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
+                num_surf(j,i)         = num_surf(j,i) + 1
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDDO
+
 
  END SUBROUTINE surface_data_output_rrd_local_mpi
@@ -4673,7 +4782,11 @@
        WRITE ( 14 )  average_count_surf
 
+       CALL wrd_write_string( 'time_dosurf_av' )
+       WRITE ( 14 )  time_dosurf_av
+
     ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
 
       CALL wrd_mpi_io( 'average_count_surf', average_count_surf )
+      CALL wrd_mpi_io( 'time_dosurf_av', time_dosurf_av )
 
     ENDIF
@@ -4691,4 +4804,24 @@
     IMPLICIT NONE
 
+    CHARACTER(LEN=3) ::  dum !< dummy string to create output-variable name
+
+    INTEGER(iwp) ::  i                      !< grid index in x-direction
+    INTEGER(iwp) ::  j                      !< grid index in y-direction
+    INTEGER(iwp) ::  l                      !< running index surface orientation
+    INTEGER(iwp) ::  m                      !< running index surface elements
+    INTEGER(iwp) ::  n                      !< counting variable
+    INTEGER(iwp) ::  nv                     !< running index over number of variables
+    INTEGER(iwp) ::  start_index_aggregated !< sum of start-index at (j,i) over all surface types
+
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  end_index           !< end index of surface data at (j,i)
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  global_start_index  !< index array for surface data (MPI-IO)
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  num_surf            !< number of surface data at (j,i)
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  start_index         !< start index of surface data at (j,i)
+
+    LOGICAL ::  surface_data_to_write  !< switch for MPI-I/O if PE has surface data to write
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surf_out  !< surface data in expected restart format
+
+
     IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
 
@@ -4699,6 +4832,134 @@
 
     ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
-
-       IF ( ALLOCATED( surfaces%var_av ) )  CALL wrd_mpi_io( 'surfaces%var_av', surfaces%var_av )
+!
+!--    Note, surface data which is written to file is organized in a different way than
+!--    the output surface data. The output surface data is a concatenated array of the 
+!--    different surface types and orientations, while the mpi-io expects surface data that
+!--    is consecutive in terms of start- and end-index, i.e. organized along the (j,i)
+!--    grid index. Hence, data need to be tranformed before it can be written to file.
+       IF ( ALLOCATED( surfaces%var_av ) )  THEN
+          ALLOCATE( end_index(nys:nyn,nxl:nxr)          )
+          ALLOCATE( num_surf(nys:nyn,nxl:nxr)           )
+          ALLOCATE( start_index(nys:nyn,nxl:nxr)        )
+          ALLOCATE( global_start_index(nys:nyn,nxl:nxr) )
+          ALLOCATE( surf_out(1:surfaces%ns)             )
+!
+!--       Determine the start and end index at each (j,i)-pair and resort the surface data
+          start_index            = 1
+          end_index              = 0
+          start_index_aggregated = 1
+          num_surf = 0
+          DO  l = 0, 1
+             DO  m = 1, surf_def_h(l)%ns
+                i = surf_def_h(l)%i(m)
+                j = surf_def_h(l)%j(m)
+                num_surf(j,i) = num_surf(j,i) + 1
+             ENDDO
+          ENDDO
+          DO  m = 1, surf_lsm_h%ns
+             i = surf_lsm_h%i(m)
+             j = surf_lsm_h%j(m)
+             num_surf(j,i) = num_surf(j,i) + 1
+          ENDDO
+          DO  m = 1, surf_usm_h%ns
+             i = surf_usm_h%i(m)
+             j = surf_usm_h%j(m)
+             num_surf(j,i) = num_surf(j,i) + 1
+          ENDDO
+
+          DO  l = 0, 3
+             DO  m = 1, surf_def_v(l)%ns
+                i = surf_def_v(l)%i(m)
+                j = surf_def_v(l)%j(m)
+                num_surf(j,i) = num_surf(j,i) + 1
+             ENDDO
+             DO  m = 1, surf_lsm_v(l)%ns
+                i = surf_lsm_v(l)%i(m)
+                j = surf_lsm_v(l)%j(m)
+                num_surf(j,i) = num_surf(j,i) + 1
+             ENDDO
+             DO  m = 1, surf_usm_v(l)%ns
+                i = surf_usm_v(l)%i(m)
+                j = surf_usm_v(l)%j(m)
+                num_surf(j,i) = num_surf(j,i) + 1
+             ENDDO
+          ENDDO
+
+          start_index = 0
+          end_index   = 0
+          start_index_aggregated = 1
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                start_index(j,i)       = start_index_aggregated
+                end_index(j,i)         = start_index(j,i) + num_surf(j,i) - 1
+                start_index_aggregated = start_index_aggregated + num_surf(j,i)
+             ENDDO
+          ENDDO
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, surface_data_to_write,         &
+                                            global_start_index )
+          CALL wrd_mpi_io( 'surfaces%start_index',        start_index        )
+          CALL wrd_mpi_io( 'surfaces%end_index',          end_index          )
+          CALL wrd_mpi_io( 'surfaces%global_start_index', global_start_index )
+
+          DO  nv = 1, dosurf_no(1)
+             n = 0
+             num_surf = 0
+             DO  l = 0, 1
+                DO  m = 1, surf_def_h(l)%ns
+                   i = surf_def_h(l)%i(m)
+                   j = surf_def_h(l)%j(m)
+                   n                                          = n + 1
+                   surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
+                   num_surf(j,i)                              = num_surf(j,i) + 1
+                ENDDO
+             ENDDO
+             DO  m = 1, surf_lsm_h%ns
+                i = surf_lsm_h%i(m)
+                j = surf_lsm_h%j(m)
+                n                                          = n + 1
+                surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
+                num_surf(j,i)                              = num_surf(j,i) + 1
+             ENDDO
+             DO  m = 1, surf_usm_h%ns
+                i = surf_usm_h%i(m)
+                j = surf_usm_h%j(m)
+                n                                          = n + 1
+                surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
+                num_surf(j,i)                              = num_surf(j,i) + 1
+             ENDDO
+
+             DO  l = 0, 3
+                DO  m = 1, surf_def_v(l)%ns
+                   i = surf_def_v(l)%i(m)
+                   j = surf_def_v(l)%j(m)
+                   n                                          = n + 1
+                   surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
+                   num_surf(j,i)                              = num_surf(j,i) + 1
+                ENDDO
+                DO  m = 1, surf_lsm_v(l)%ns
+                   i = surf_lsm_v(l)%i(m)
+                   j = surf_lsm_v(l)%j(m)
+                   n                                          = n + 1
+                   surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
+                   num_surf(j,i)                              = num_surf(j,i) + 1
+                ENDDO
+                DO  m = 1, surf_usm_v(l)%ns
+                   i = surf_usm_v(l)%i(m)
+                   j = surf_usm_v(l)%j(m)
+                   n                                          = n + 1
+                   surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
+                   num_surf(j,i)                              = num_surf(j,i) + 1
+                ENDDO
+             ENDDO
+
+             IF ( nv < 10                     )  WRITE( dum, '(I1)') nv
+             IF ( nv < 100   .AND.  nv >= 10  )  WRITE( dum, '(I2)') nv
+             IF ( nv < 1000  .AND.  nv >= 100 )  WRITE( dum, '(I3)') nv
+
+             CALL wrd_mpi_io_surface( 'surfaces%var_av' // TRIM( dum ), surf_out )
+          ENDDO
+
+       ENDIF
 
     ENDIF
