Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 4513)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 4514)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! Enable output of qsurf and ssurf
+! 
+! 4509 2020-04-26 15:57:55Z raasch
 ! file re-formatted to follow the PALM coding standard
 ! 
@@ -284,4 +287,14 @@
              ENDIF
 
+         CASE ( 'qsurf*' )
+             IF ( ALLOCATED( qsurf_av ) )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qsurf_av(j,i) = qsurf_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+                CALL exchange_horiz_2d( qsurf_av )
+             ENDIF
+
          CASE ( 'r_a*' )
              IF ( ALLOCATED( r_a_av ) )  THEN
@@ -315,4 +328,14 @@
              ENDIF
 
+          CASE ( 'ssurf*' )
+             IF ( ALLOCATED( ssurf_av ) )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      ssurf_av(j,i) = ssurf_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+                CALL exchange_horiz_2d( ssurf_av )
+             ENDIF
+
           CASE ( 'ssws*' )
              IF ( ALLOCATED( ssws_av ) )  THEN
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 4513)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 4514)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enable output of qsurf and ssurf
+! 
+! 4513 2020-04-30 13:45:47Z raasch
 ! unused modules removed
 ! 
@@ -2478,7 +2481,7 @@
              CONTINUE
 
-          CASE ( 'ghf*', 'lwp*', 'ol*', 'qsws*', 'r_a*',                       &
-                 'shf*', 'ssws*', 't*', 'tsurf*', 'us*',                       &
-                 'z0*', 'z0h*', 'z0q*' )
+          CASE ( 'ghf*', 'lwp*', 'ol*', 'qsurf*', 'qsws*', 'r_a*',             &
+                 'shf*', 'ssurf*', 'ssws*', 't*', 'tsurf*', 'us*',             &
+                 'z0*', 'z0h*', 'z0q*'  )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                 message_string = 'illegal value for data_output: "' //         &
@@ -2523,7 +2526,9 @@
              IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
              IF ( TRIM( var ) == 'ol*'    )  unit = 'm'
+             IF ( TRIM( var ) == 'qsurf*' )  unit = 'kg/kg'
              IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
              IF ( TRIM( var ) == 'r_a*'   )  unit = 's/m'
              IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
+             IF ( TRIM( var ) == 'ssurf*' )  unit = 'kg/kg'
              IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'
              IF ( TRIM( var ) == 't*'     )  unit = 'K'
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 4513)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 4514)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enable output of qsurf and ssurf
+! 
+! 4500 2020-04-17 10:12:45Z suehring
 ! Unify output conversion of sensible and latent heat flux
 ! 
@@ -123,5 +126,7 @@
     USE indices,                                                               &
         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,       &
-               nzb, nzt, wall_flags_total_0
+               nzb, nzt,                                                       &
+               topo_top_ind,                                                   &
+               wall_flags_total_0
 
     USE kinds
@@ -701,4 +706,39 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'qsurf*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  m = 1, surf_def_h(0)%ns
+                      i                   = surf_def_h(0)%i(m)
+                      j                   = surf_def_h(0)%j(m)
+                      local_pf(i,j,nzb+1) = surf_def_h(0)%q_surface(m)
+                   ENDDO
+
+                   DO  m = 1, surf_lsm_h%ns
+                      i                   = surf_lsm_h%i(m) 
+                      j                   = surf_lsm_h%j(m)
+                      local_pf(i,j,nzb+1) = surf_lsm_h%q_surface(m)
+                   ENDDO
+
+                   DO  m = 1, surf_usm_h%ns
+                      i                   = surf_usm_h%i(m)
+                      j                   = surf_usm_h%j(m)
+                      local_pf(i,j,nzb+1) = surf_usm_h%q_surface(m)
+                   ENDDO
+
+                ELSE
+                   IF ( .NOT. ALLOCATED( qsurf_av ) ) THEN
+                      ALLOCATE( qsurf_av(nysg:nyng,nxlg:nxrg) )
+                      qsurf_av = REAL( fill_value, KIND = wp )
+                   ENDIF
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         local_pf(i,j,nzb+1) = qsurf_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted       = .TRUE.
+                two_d          = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'qsws*_xy' )        ! 2d-array
                 IF ( av == 0 ) THEN
@@ -845,5 +885,28 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
-                
+
+             CASE ( 'ssurf*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         k = topo_top_ind(j,i,0)
+                         local_pf(i,j,nzb+1) = s(k+k,j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   IF ( .NOT. ALLOCATED( ssurf_av ) ) THEN
+                      ALLOCATE( ssurf_av(nysg:nyng,nxlg:nxrg) )
+                      ssurf_av = REAL( fill_value, KIND = wp )
+                   ENDIF
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         local_pf(i,j,nzb+1) = ssurf_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted       = .TRUE.
+                two_d          = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'ssws*_xy' )        ! 2d-array
                 IF ( av == 0 ) THEN
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 4513)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 4514)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Add possibility to initialize surface sensible and latent heat fluxes via
+! a static driver.
+! 
+! 4493 2020-04-10 09:49:43Z pavelkrc
 ! Overwrite u_init, v_init, pt_init, q_init and s_init with hom for all
 ! cyclic_fill-cases, not only for turbulent_inflow = .TRUE.
@@ -1599,5 +1603,5 @@
 !--       Initialize roughness length. Note, z0 will be only initialized at    
 !--       default-type surfaces. At natural or urban z0 is implicitly          
-!--       initialized bythe respective parameter lists.                        
+!--       initialized by the respective parameter lists.                        
 !--       Initialize horizontal surface elements.                              
           CALL init_single_surface_properties( surf_def_h(0)%z0,               &
@@ -1624,7 +1628,83 @@
        ENDIF
 !
-!--    Additional variables, e.g. shf, qsws, etc, can be initialized the 
+!--    Input surface sensible heat flux. 
+       IF ( check_existence( vars_pids, 'shf' ) )  THEN
+!
+!--       Read _FillValue attribute
+          CALL get_attribute( pids_id, char_fill, tmp_2d%fill,                 &
+                              .FALSE., 'shf' )
+!
+!--       Read variable
+          CALL get_variable( pids_id, 'shf', tmp_2d%var,                       &
+                             nxl, nxr, nys, nyn )
+!
+!--       Initialize heat flux. Note, shf will be only initialized at
+!--       default-type surfaces. At natural or urban shf is implicitly
+!--       initialized by the respective parameter lists.
+!--       Initialize horizontal surface elements.
+          CALL init_single_surface_properties( surf_def_h(0)%shf,              &
+                                               tmp_2d%var,                     &
+                                               surf_def_h(0)%ns,               &
+                                               tmp_2d%fill,                    &
+                                               surf_def_h(0)%i,                &
+                                               surf_def_h(0)%j )
+!
+!--       Initialize heat flux also at vertical surface elements.
+!--       Note, the actual 2D input arrays are only defined on the
+!--       subdomain. Therefore, pass the index arrays with their respective
+!--       offset values.
+          DO  l = 0, 3
+             CALL init_single_surface_properties(                              &
+                                         surf_def_v(l)%shf,                    &
+                                         tmp_2d%var,                           &
+                                         surf_def_v(l)%ns,                     &
+                                         tmp_2d%fill,                          &
+                                         surf_def_v(l)%i + surf_def_v(l)%ioff, &
+                                         surf_def_v(l)%j + surf_def_v(l)%joff )
+          ENDDO
+
+       ENDIF
+!
+!--    Input surface sensible heat flux. 
+       IF ( check_existence( vars_pids, 'qsws' ) )  THEN
+!
+!--       Read _FillValue attribute
+          CALL get_attribute( pids_id, char_fill, tmp_2d%fill,                 &
+                              .FALSE., 'qsws' )
+!
+!--       Read variable
+          CALL get_variable( pids_id, 'qsws', tmp_2d%var,                      &
+                             nxl, nxr, nys, nyn )
+!
+!--       Initialize latent heat flux. Note, qsws will be only initialized at
+!--       default-type surfaces. At natural or urban qsws is implicitly
+!--       initialized by the respective parameter lists.
+!--       Initialize horizontal surface elements.
+          CALL init_single_surface_properties( surf_def_h(0)%qsws,             &
+                                               tmp_2d%var,                     &
+                                               surf_def_h(0)%ns,               &
+                                               tmp_2d%fill,                    &
+                                               surf_def_h(0)%i,                &
+                                               surf_def_h(0)%j )
+!
+!--       Initialize latent heat flux also at vertical surface elements.
+!--       Note, the actual 2D input arrays are only defined on the
+!--       subdomain. Therefore, pass the index arrays with their respective
+!--       offset values.
+          DO  l = 0, 3
+             CALL init_single_surface_properties(                              &
+                                         surf_def_v(l)%qsws,                   &
+                                         tmp_2d%var,                           &
+                                         surf_def_v(l)%ns,                     &
+                                         tmp_2d%fill,                          &
+                                         surf_def_v(l)%i + surf_def_v(l)%ioff, &
+                                         surf_def_v(l)%j + surf_def_v(l)%joff )
+          ENDDO
+
+       ENDIF
+!
+!--    Additional variables, can be initialized the 
 !--    same way.
-                           
+
 !
 !--    Finally, close the input file and deallocate temporary arrays
Index: palm/trunk/SOURCE/module_interface.f90
===================================================================
--- palm/trunk/SOURCE/module_interface.f90	(revision 4513)
+++ palm/trunk/SOURCE/module_interface.f90	(revision 4514)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added global restart routines for plant-canopy model
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
 ! 
@@ -442,5 +445,7 @@
                pcm_3d_data_averaging,                                          &
                pcm_data_output_3d,                                             &
+               pcm_rrd_global,                                                 &
                pcm_rrd_local,                                                  &
+               pcm_wrd_global,                                                 &
                pcm_wrd_local
 
@@ -1722,4 +1727,5 @@
     IF ( .NOT. found )  CALL gust_rrd_global( found ) ! ToDo: change interface to pass variable
     IF ( .NOT. found )  CALL lpm_rrd_global( found ) ! ToDo: change interface to pass variable
+    IF ( .NOT. found )  CALL pcm_rrd_global( found )
     IF ( .NOT. found )  CALL ocean_rrd_global( found ) ! ToDo: change interface to pass variable
     IF ( .NOT. found )  CALL stg_rrd_global ( found ) ! ToDo: change interface to pass variable
@@ -1752,4 +1758,5 @@
     IF ( gust_module_enabled )  CALL gust_rrd_global
     IF ( particle_advection )   CALL lpm_rrd_global
+    IF ( plant_canopy )         CALL pcm_rrd_global
     IF ( ocean_mode )           CALL ocean_rrd_global
     IF ( syn_turb_gen )         CALL stg_rrd_global
@@ -1782,4 +1789,5 @@
     IF ( gust_module_enabled )  CALL gust_wrd_global
     IF ( particle_advection )   CALL lpm_wrd_global
+    IF ( plant_canopy )         CALL pcm_wrd_global
     IF ( ocean_mode )           CALL ocean_wrd_global
     IF ( syn_turb_gen )         CALL stg_wrd_global
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 4513)
+++ palm/trunk/SOURCE/modules.f90	(revision 4514)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! +qsurf_av, ssurf_av
+! 
+! 4505 2020-04-20 15:37:15Z schwenkel
 ! Add flag for saturation check
 ! 
@@ -471,6 +474,8 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av                 !< avg. liquid water path
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ol_av                  !< avg. Obukhov length
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qsurf_av               !< avg. surface mixing ratio
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qsws_av                !< avg. surface moisture flux
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  r_a_av                 !< avg. resistance
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ssurf_av               !< avg. surface passive scalar
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ssws_av                !< avg. surface scalar flux
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  shf_av                 !< avg. surface heat flux
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4513)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4514)
@@ -27,4 +27,11 @@
 ! -----------------
 ! $Id$
+! - Bugfix in output of pcm_heatrate_av in a restart run. In order to fix this, 
+!   pch_index is now output for a restart run. Therefore, define global restart
+!   routines.
+! - Error message number renamed and check for PA0505 revised in order to also
+!   consider natural surfaces with plant-canopy.
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
 !
@@ -168,4 +175,5 @@
 ! @todo - precalculate constant terms in pcm_calc_transpiration_rate
 ! @todo - unify variable names (pcm_, pc_, ...)
+! @todo - get rid-off dependency on radiation model
 !------------------------------------------------------------------------------!
  MODULE plant_canopy_model_mod
@@ -187,4 +195,5 @@
               dz,                                                              &
               humidity,                                                        &
+              land_surface,                                                    &
               length,                                                          &
               message_string,                                                  &
@@ -225,5 +234,6 @@
 
     USE restart_data_mpi_io_mod,                                                                   &
-        ONLY:  wrd_mpi_io
+        ONLY:  rrd_mpi_io,                                                                         &
+               wrd_mpi_io
 
     USE surface_mod,                                                           &
@@ -294,6 +304,8 @@
            pcm_init,                                                          &
            pcm_parin,                                                         &
+           pcm_rrd_global,                                                    &
            pcm_rrd_local,                                                     &
            pcm_tendency,                                                      &
+           pcm_wrd_global,                                                    &
            pcm_wrd_local
 
@@ -349,4 +361,9 @@
     END INTERFACE pcm_rrd_local
 
+    INTERFACE pcm_rrd_global
+       MODULE PROCEDURE pcm_rrd_global_ftn
+       MODULE PROCEDURE pcm_rrd_global_mpi
+    END INTERFACE pcm_rrd_global
+
     INTERFACE pcm_tendency
        MODULE PROCEDURE pcm_tendency
@@ -357,4 +374,8 @@
        MODULE PROCEDURE pcm_wrd_local
     END INTERFACE pcm_wrd_local
+
+    INTERFACE pcm_wrd_global
+       MODULE PROCEDURE pcm_wrd_global
+    END INTERFACE pcm_wrd_global
 
 
@@ -483,8 +504,17 @@
 
        CASE ( 'pcm_heatrate' )
-          IF ( cthf == 0.0_wp  .AND. .NOT.  urban_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+!
+!--       Output of heatrate can be only done if it is explicitely set by cthf,
+!--       or parametrized by absorption of radiation. The latter, however, is
+!--       only available if radiation_interactions are on. Note, these are 
+!--       enabled if land-surface or urban-surface is switched-on. Using 
+!--       radiation_interactions_on directly is not possible since it belongs
+!--       to the radition_model, which in turn depends on the plant-canopy model,
+!--       creating circular dependencies.
+          IF ( cthf == 0.0_wp  .AND. (  .NOT.  urban_surface  .AND.             &
+                                        .NOT.  land_surface ) )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //      &
                               'res setting of parameter cthf /= 0.0'
-             CALL message( 'pcm_check_data_output', 'PA1000', 1, 2, 0, 6, 0 )
+             CALL message( 'pcm_check_data_output', 'PA0505', 1, 2, 0, 6, 0 )
           ENDIF
           unit = 'K s-1'
@@ -1620,4 +1650,39 @@
 ! Description:
 ! ------------
+!> Read module-specific global restart data (Fortran binary format).
+!------------------------------------------------------------------------------!
+    SUBROUTINE pcm_rrd_global_ftn( found )
+
+       LOGICAL, INTENT(OUT)  ::  found
+
+       found = .TRUE.
+
+       SELECT CASE ( restart_string(1:length) )
+
+          CASE ( 'pch_index' )
+             READ ( 13 )  pch_index
+
+          CASE DEFAULT
+
+             found = .FALSE.
+
+       END SELECT
+
+    END SUBROUTINE pcm_rrd_global_ftn
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read module-specific global restart data (MPI-IO).
+!------------------------------------------------------------------------------!
+    SUBROUTINE pcm_rrd_global_mpi
+
+       CALL rrd_mpi_io( 'pch_index', pch_index )
+
+    END SUBROUTINE pcm_rrd_global_mpi
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Subroutine reads local (subdomain) restart data
 !------------------------------------------------------------------------------!
@@ -2437,4 +2502,24 @@
 ! Description:
 ! ------------
+!> Subroutine writes global restart data
+!------------------------------------------------------------------------------!
+    SUBROUTINE pcm_wrd_global
+
+       IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
+
+          CALL wrd_write_string( 'pch_index' )
+          WRITE ( 14 )  pch_index
+
+       ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' )  THEN
+
+          CALL wrd_mpi_io( 'pch_index', pch_index )
+
+       ENDIF
+
+    END SUBROUTINE pcm_wrd_global
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Subroutine writes local (subdomain) restart data
 !------------------------------------------------------------------------------!
Index: palm/trunk/SOURCE/read_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 4513)
+++ palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 4514)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! qsurf and ssurf added
+! 
+! 4498 2020-04-15 14:26:31Z raasch
 ! argument removed from rd_mpi_io_open
 ! 
@@ -1797,4 +1800,12 @@
                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+                CASE ( 'qsurf_av' )
+                   IF ( .NOT. ALLOCATED( qsurf_av ) )  THEN
+                      ALLOCATE( qsurf_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsurf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                             &
+                      tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'qsws_av' )
                    IF ( .NOT. ALLOCATED( qsws_av ) )  THEN
@@ -1859,4 +1870,12 @@
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
                    shf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =                              &
+                      tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'ssurf_av' )
+                   IF ( .NOT. ALLOCATED( ssurf_av ) )  THEN
+                      ALLOCATE( ssurf_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   ssurf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                             &
                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4513)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4514)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enable output of qsurf and ssurf
+! 
+! 4442 2020-03-04 19:21:13Z suehring
 ! Change order of dimension in surface array %frac to allow for better 
 ! vectorization.
@@ -77,6 +80,7 @@
     USE averaging,                                                             &
         ONLY:  e_av, ghf_av, lpt_av, lwp_av, ol_av, p_av, pc_av, pr_av, pt_av, &
-               q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsws_av,               &
-               qv_av, r_a_av, s_av, shf_av, ssws_av, ts_av, tsurf_av, u_av,    &
+               q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsurf_av, qsws_av,     &
+               qv_av, r_a_av, s_av, shf_av, ssurf_av,                          &
+               ssws_av, ts_av, tsurf_av, u_av,                                 &
                us_av, v_av, vpt_av, w_av, z0_av, z0h_av, z0q_av
 
@@ -95,5 +99,6 @@
 
     USE indices,                                                               &
-        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt 
+        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,           &
+               topo_top_ind
 
     USE kinds
@@ -232,4 +237,10 @@
                 ql_vp_av = 0.0_wp
 
+             CASE ( 'qsurf*' )
+                IF ( .NOT. ALLOCATED( qsurf_av ) )  THEN
+                   ALLOCATE( qsurf_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsurf_av = 0.0_wp
+
              CASE ( 'qsws*' )
                 IF ( .NOT. ALLOCATED( qsws_av ) )  THEN
@@ -261,5 +272,11 @@
                 ENDIF
                 shf_av = 0.0_wp
-                
+
+             CASE ( 'ssurf*' )
+                IF ( .NOT. ALLOCATED( ssurf_av ) )  THEN
+                   ALLOCATE( ssurf_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                ssurf_av = 0.0_wp
+
              CASE ( 'ssws*' )
                 IF ( .NOT. ALLOCATED( ssws_av ) )  THEN
@@ -587,4 +604,32 @@
              ENDIF
 
+          CASE ( 'qsurf*' )
+             IF ( ALLOCATED( qsurf_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      match_def = surf_def_h(0)%start_index(j,i) <=            &
+                                  surf_def_h(0)%end_index(j,i)
+                      match_lsm = surf_lsm_h%start_index(j,i) <=               &
+                                  surf_lsm_h%end_index(j,i)
+                      match_usm = surf_usm_h%start_index(j,i) <=               &
+                                  surf_usm_h%end_index(j,i)
+
+                      IF ( match_def )  THEN
+                         m = surf_def_h(0)%end_index(j,i)
+                         qsurf_av(j,i) = qsurf_av(j,i) +                       &
+                                         surf_def_h(0)%q_surface(m)
+                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
+                         m = surf_lsm_h%end_index(j,i)
+                         qsurf_av(j,i) = qsurf_av(j,i) +                       &
+                                         surf_lsm_h%q_surface(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         qsurf_av(j,i) = qsurf_av(j,i) +                       &
+                                         surf_usm_h%q_surface(m)
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDIF
+
           CASE ( 'qsws*' )
 !
@@ -703,4 +748,14 @@
              ENDIF
 
+          CASE ( 'ssurf*' )
+             IF ( ALLOCATED( ssurf_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      k = topo_top_ind(j,i,0)
+                      ssurf_av(j,i) = ssurf_av(j,i) + s(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
+
           CASE ( 'ssws*' )
              IF ( ALLOCATED( ssws_av ) ) THEN 
Index: palm/trunk/SOURCE/write_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 4513)
+++ palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 4514)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! qsurf and ssurf added
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
 ! 
@@ -1254,4 +1257,9 @@
           ENDIF
 
+          IF ( ALLOCATED( qsurf_av ) )  THEN
+             CALL wrd_write_string( 'qsurf_av' )
+             WRITE ( 14 )  qsurf_av
+          ENDIF
+
           IF ( ALLOCATED( qsws_av ) )  THEN
              CALL wrd_write_string( 'qsws_av' )
@@ -1269,4 +1277,9 @@
              CALL wrd_write_string( 's_av' )
              WRITE ( 14 )  s_av
+          ENDIF
+
+          IF ( ALLOCATED( ssurf_av ) )  THEN
+             CALL wrd_write_string( 'ssurf_av' )
+             WRITE ( 14 )  ssurf_av
           ENDIF
 
@@ -1491,4 +1504,5 @@
           ENDIF
 
+          IF ( ALLOCATED( qsurf_av ) )  CALL wrd_mpi_io( 'qsurf_av', qsurf_av )
           IF ( ALLOCATED( qsws_av ) )  CALL wrd_mpi_io( 'qsws_av', qsws_av )
 
@@ -1499,4 +1513,5 @@
           CALL wrd_mpi_io( 's', s )
           IF ( ALLOCATED( s_av ) )  CALL wrd_mpi_io( 's_av',  s_av )
+          IF ( ALLOCATED( ssurf_av ) )  CALL wrd_mpi_io( 'ssurf_av', ssurf_av )
           IF ( ALLOCATED( ssws_av ) )  CALL wrd_mpi_io( 'ssws_av', ssws_av )
 
