Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4901)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1)
@@ -1,72 +1,12 @@
-!> @file sum_up_3d_data.f90
-!--------------------------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
-! Public License as published by the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
-! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
-! Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with PALM. If not, see
-! <http://www.gnu.org/licenses/>.
-!
-! Copyright 1997-2021 Leibniz Universitaet Hannover
-!--------------------------------------------------------------------------------------------------!
-!
-!
-! Current revisions:
+ SUBROUTINE sum_up_3d_data
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
 ! -----------------
-!
+! 
 !
 ! Former revisions:
 ! -----------------
-! $Id$
-! Implementation of downward facing USM and LSM surfaces
-!
-! 4591 2020-07-06 15:56:08Z raasch
-! File re-formatted to follow the PALM coding standard
-!
-! 4516 2020-04-30 16:55:10Z suehring
-! Remove double index
-!
-! 4514 2020-04-30 16:29:59Z suehring
-! 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.
-!
-! 4441 2020-03-04 19:20:35Z suehring
-! Move 2-m potential temperature output to diagnostic_output_quantities
-!
-! 4182 2019-08-22 15:20:23Z scharf
-! Corrected "Former revisions" section
-!
-! 4048 2019-06-21 21:00:21Z knoop
-! Moved tcm_3d_data_averaging to module_interface
-!
-! 4039 2019-06-18 10:32:41Z suehring
-! Modularize diagnostic output
-!
-! 3994 2019-05-22 18:08:09Z suehring
-! Output of turbulence intensity added
-!
-! 3943 2019-05-02 09:50:41Z maronga
-! Added output of qsws_av for green roofs.
-!
-! 3933 2019-04-25 12:33:20Z kanani
-! Formatting
-!
-! 3773 2019-03-01 08:56:57Z maronga
-! Added output of theta_2m*_xy_av
-!
-! 3761 2019-02-25 15:31:42Z raasch
-! unused variables removed
-!
-! 3655 2019-01-07 16:51:22Z knoop
-! Implementation of the PALM module interface
-!
+! $Log: sum_up_3d_data.f90,v $
 ! Revision 1.1  2006/02/23 12:55:23  raasch
 ! Initial revision
@@ -75,136 +15,22 @@
 ! Description:
 ! ------------
-!> Sum-up the values of 3d-arrays. The real averaging is later done in routine average_3d_data.
-!--------------------------------------------------------------------------------------------------!
- SUBROUTINE sum_up_3d_data
-
-
-    USE arrays_3d,                                                                                 &
-        ONLY:  dzw,                                                                                &
-               d_exner,                                                                            &
-               e,                                                                                  &
-               heatflux_output_conversion,                                                         &
-               p,                                                                                  &
-               pt,                                                                                 &
-               q,                                                                                  &
-               ql,                                                                                 &
-               ql_c,                                                                               &
-               ql_v,                                                                               &
-               s,                                                                                  &
-               u,                                                                                  &
-               v,                                                                                  &
-               vpt,                                                                                &
-               w,                                                                                  &
-               waterflux_output_conversion
-
-    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,                                                                           &
-               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
-
-    USE basic_constants_and_equations_mod,                                                         &
-        ONLY:  c_p,                                                                                &
-               lv_d_cp,                                                                            &
-               l_v
-
-    USE bulk_cloud_model_mod,                                                                      &
-        ONLY:  bulk_cloud_model
-
-    USE control_parameters,                                                                        &
-        ONLY:  average_count_3d,                                                                   &
-               doav,                                                                               &
-               doav_n,                                                                             &
-               rho_surface,                                                                        &
-               urban_surface,                                                                      &
-               varnamelength
-
-    USE cpulog,                                                                                    &
-        ONLY:  cpu_log,                                                                            &
-               log_point
-
-    USE indices,                                                                                   &
-        ONLY:  nxl,                                                                                &
-               nxlg,                                                                               &
-               nxr,                                                                                &
-               nxrg,                                                                               &
-               nyn,                                                                                &
-               nyng,                                                                               &
-               nys,                                                                                &
-               nysg,                                                                               &
-               nzb,                                                                                &
-               nzt,                                                                                &
-               topo_top_ind
-
-    USE kinds
-
-    USE module_interface,                                                                          &
-        ONLY:  module_interface_3d_data_averaging
-
-    USE particle_attributes,                                                                       &
-        ONLY:  grid_particles,                                                                     &
-               number_of_particles,                                                                &
-               particles,                                                                          &
-               prt_count
-
-    USE surface_mod,                                                                               &
-        ONLY:  ind_pav_green,                                                                      &
-               ind_veg_wall,                                                                       &
-               ind_wat_win,                                                                        &
-               surf_def_h,                                                                         &
-               surf_lsm_h,                                                                         &
-               surf_usm_h
-
-    USE urban_surface_mod,                                                                         &
-        ONLY:  usm_3d_data_averaging
-
+! Sum-up the values of 3d-arrays. The real averaging is later done in routine
+! average_3d_data. 
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    USE averaging
+    USE cloud_parameters
+    USE control_parameters
+    USE cpulog
+    USE indices
+    USE interfaces
+    USE particle_attributes
 
     IMPLICIT NONE
 
-    CHARACTER(LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
-
-    INTEGER(iwp) ::  i   !< grid index x direction
-    INTEGER(iwp) ::  ii  !< running index
-    INTEGER(iwp) ::  j   !< grid index y direction
-    INTEGER(iwp) ::  k   !< grid index x direction
-    INTEGER(iwp) ::  m   !< running index over surfacle elements
-    INTEGER(iwp) ::  n   !< running index over number of particles per grid box
-
-    LOGICAL ::  match_def  !< flag indicating default-type surface
-    LOGICAL ::  match_lsm  !< flag indicating natural-type surface
-    LOGICAL ::  match_usm  !< flag indicating urban-type surface
-
-    REAL(wp) ::  mean_r  !< mean-particle radius witin grid box
-    REAL(wp) ::  s_r2    !< mean-particle radius witin grid box to the power of two
-    REAL(wp) ::  s_r3    !< mean-particle radius witin grid box to the power of three
-
+    INTEGER ::  i, ii, j, k, n, psi
+
+    REAL    ::  mean_r, s_r3, s_r4
 
 
@@ -212,214 +38,132 @@
 
 !
-!-- Allocate and initialize the summation arrays if called for the very first time or the first time
-!-- after average_3d_data has been called (some or all of the arrays may have been already allocated
-!-- in rrd_local)
+!-- Allocate and initialize the summation arrays if called for the very first
+!-- time or the first time after average_3d_data has been called 
+!-- (some or all of the arrays may have been already allocated
+!-- in read_3d_binary)
     IF ( average_count_3d == 0 )  THEN
 
        DO  ii = 1, doav_n
 
-          trimvar = TRIM( doav(ii) )
-
-          SELECT CASE ( trimvar )
-
-             CASE ( 'ghf*' )
-                IF ( .NOT. ALLOCATED( ghf_av ) )  THEN
-                   ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                ghf_av = 0.0_wp
+          SELECT CASE ( TRIM( doav(ii) ) )
 
              CASE ( 'e' )
                 IF ( .NOT. ALLOCATED( e_av ) )  THEN
-                   ALLOCATE( e_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                e_av = 0.0_wp
-
-             CASE ( 'thetal' )
-                IF ( .NOT. ALLOCATED( lpt_av ) )  THEN
-                   ALLOCATE( lpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                lpt_av = 0.0_wp
+                   ALLOCATE( e_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                e_av = 0.0
 
              CASE ( 'lwp*' )
                 IF ( .NOT. ALLOCATED( lwp_av ) )  THEN
-                   ALLOCATE( lwp_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                lwp_av = 0.0_wp
-
-             CASE ( 'ol*' )
-                IF ( .NOT. ALLOCATED( ol_av ) )  THEN
-                   ALLOCATE( ol_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                ol_av = 0.0_wp
+                   ALLOCATE( lwp_av(nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                lwp_av = 0.0
 
              CASE ( 'p' )
                 IF ( .NOT. ALLOCATED( p_av ) )  THEN
-                   ALLOCATE( p_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                p_av = 0.0_wp
+                   ALLOCATE( p_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                p_av = 0.0
 
              CASE ( 'pc' )
                 IF ( .NOT. ALLOCATED( pc_av ) )  THEN
-                   ALLOCATE( pc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                pc_av = 0.0_wp
+                   ALLOCATE( pc_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                pc_av = 0.0
 
              CASE ( 'pr' )
                 IF ( .NOT. ALLOCATED( pr_av ) )  THEN
-                   ALLOCATE( pr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                pr_av = 0.0_wp
-
-             CASE ( 'theta' )
+                   ALLOCATE( pr_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                pr_av = 0.0
+
+             CASE ( 'pt' )
                 IF ( .NOT. ALLOCATED( pt_av ) )  THEN
-                   ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                pt_av = 0.0_wp
+                   ALLOCATE( pt_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                pt_av = 0.0
 
              CASE ( 'q' )
                 IF ( .NOT. ALLOCATED( q_av ) )  THEN
-                   ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                q_av = 0.0_wp
+                   ALLOCATE( q_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                q_av = 0.0
 
              CASE ( 'ql' )
                 IF ( .NOT. ALLOCATED( ql_av ) )  THEN
-                   ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                ql_av = 0.0_wp
+                   ALLOCATE( ql_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                ql_av = 0.0
 
              CASE ( 'ql_c' )
                 IF ( .NOT. ALLOCATED( ql_c_av ) )  THEN
-                   ALLOCATE( ql_c_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                ql_c_av = 0.0_wp
+                   ALLOCATE( ql_c_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                ql_c_av = 0.0
 
              CASE ( 'ql_v' )
                 IF ( .NOT. ALLOCATED( ql_v_av ) )  THEN
-                   ALLOCATE( ql_v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                ql_v_av = 0.0_wp
+                   ALLOCATE( ql_v_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                ql_v_av = 0.0
 
              CASE ( 'ql_vp' )
                 IF ( .NOT. ALLOCATED( ql_vp_av ) )  THEN
-                   ALLOCATE( ql_vp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                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
-                   ALLOCATE( qsws_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                qsws_av = 0.0_wp
+                   ALLOCATE( ql_vp_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                ql_vp_av = 0.0
 
              CASE ( 'qv' )
                 IF ( .NOT. ALLOCATED( qv_av ) )  THEN
-                   ALLOCATE( qv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                qv_av = 0.0_wp
-
-             CASE ( 'r_a*' )
-                IF ( .NOT. ALLOCATED( r_a_av ) )  THEN
-                   ALLOCATE( r_a_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                r_a_av = 0.0_wp
+                   ALLOCATE( qv_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                qv_av = 0.0
 
              CASE ( 's' )
                 IF ( .NOT. ALLOCATED( s_av ) )  THEN
-                   ALLOCATE( s_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                s_av = 0.0_wp
-
-             CASE ( 'shf*' )
-                IF ( .NOT. ALLOCATED( shf_av ) )  THEN
-                   ALLOCATE( shf_av(nysg:nyng,nxlg:nxrg) )
-                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
-                   ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                ssws_av = 0.0_wp
+                   ALLOCATE( s_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                s_av = 0.0
 
              CASE ( 't*' )
                 IF ( .NOT. ALLOCATED( ts_av ) )  THEN
-                   ALLOCATE( ts_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                ts_av = 0.0_wp
-
-             CASE ( 'tsurf*' )
-                IF ( .NOT. ALLOCATED( tsurf_av ) )  THEN
-                   ALLOCATE( tsurf_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                tsurf_av = 0.0_wp
+                   ALLOCATE( ts_av(nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                ts_av = 0.0
 
              CASE ( 'u' )
                 IF ( .NOT. ALLOCATED( u_av ) )  THEN
-                   ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                u_av = 0.0_wp
-
-             CASE ( 'us*' )
+                   ALLOCATE( u_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                u_av = 0.0
+
+             CASE ( 'u*' )
                 IF ( .NOT. ALLOCATED( us_av ) )  THEN
-                   ALLOCATE( us_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                us_av = 0.0_wp
+                   ALLOCATE( us_av(nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                us_av = 0.0
 
              CASE ( 'v' )
                 IF ( .NOT. ALLOCATED( v_av ) )  THEN
-                   ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                v_av = 0.0_wp
-
-             CASE ( 'thetav' )
+                   ALLOCATE( v_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                v_av = 0.0
+
+             CASE ( 'vpt' )
                 IF ( .NOT. ALLOCATED( vpt_av ) )  THEN
-                   ALLOCATE( vpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                vpt_av = 0.0_wp
+                   ALLOCATE( vpt_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                vpt_av = 0.0
 
              CASE ( 'w' )
                 IF ( .NOT. ALLOCATED( w_av ) )  THEN
-                   ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                w_av = 0.0_wp
-
-             CASE ( 'z0*' )
-                IF ( .NOT. ALLOCATED( z0_av ) )  THEN
-                   ALLOCATE( z0_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                z0_av = 0.0_wp
-
-             CASE ( 'z0h*' )
-                IF ( .NOT. ALLOCATED( z0h_av ) )  THEN
-                   ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                z0h_av = 0.0_wp
-
-             CASE ( 'z0q*' )
-                IF ( .NOT. ALLOCATED( z0q_av ) )  THEN
-                   ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                z0q_av = 0.0_wp
-
+                   ALLOCATE( w_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                w_av = 0.0
 
              CASE DEFAULT
-
-!
-!--             Allocating and initializing data arrays for all other modules
-                CALL module_interface_3d_data_averaging( 'allocate', trimvar )
-
+!
+!--             User-defined quantity
+                CALL user_3d_data_averaging( 'allocate', doav(ii) )
 
           END SELECT
@@ -433,569 +177,201 @@
     DO  ii = 1, doav_n
 
-       trimvar = TRIM( doav(ii) )
 !
 !--    Store the array chosen on the temporary array.
-       SELECT CASE ( trimvar )
-
-          CASE ( 'ghf*' )
-             IF ( ALLOCATED( ghf_av ) ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-!
-!--                   Check whether grid point is a natural- or urban-type surface.
-                      match_lsm = surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-!
-!--                   In order to avoid double-counting of surface properties, always assume that
-!--                   natural-type surfaces are below urban type surfaces, e.g. in case of bridges.
-!--                   Further, take only the last suface element, i.e. the uppermost surface which
-!--                   would be visible from above
-                      IF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         ghf_av(j,i) = ghf_av(j,i) + surf_lsm_h(0)%ghf(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         ghf_av(j,i) = ghf_av(j,i) + surf_usm_h(0)%frac(m,ind_veg_wall)  *            &
-                                                     surf_usm_h(0)%wghf_eb(m)        +                &
-                                                     surf_usm_h(0)%frac(m,ind_pav_green) *            &
-                                                     surf_usm_h(0)%wghf_eb_green(m)  +                &
-                                                     surf_usm_h(0)%frac(m,ind_wat_win)   *            &
-                                                     surf_usm_h(0)%wghf_eb_window(m)
+       SELECT CASE ( TRIM( doav(ii) ) )
+
+          CASE ( 'e' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      e_av(k,j,i) = e_av(k,j,i) + e(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'lwp*' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i) * &
+                                                    dzw(1:nzt+1) )
+                ENDDO
+             ENDDO
+
+          CASE ( 'p' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      p_av(k,j,i) = p_av(k,j,i) + p(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'pc' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nzt+1
+                      pc_av(k,j,i) = pc_av(k,j,i) + prt_count(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'pr' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nzt+1
+                      psi = prt_start_index(k,j,i)
+                      s_r3 = 0.0
+                      s_r4 = 0.0
+                      DO  n = psi, psi+prt_count(k,j,i)-1
+                         s_r3 = s_r3 + particles(n)%radius**3
+                         s_r4 = s_r4 + particles(n)%radius**4
+                      ENDDO
+                      IF ( s_r3 /= 0.0 )  THEN
+                         mean_r = s_r4 / s_r3
+                      ELSE
+                         mean_r = 0.0
                       ENDIF
+                      pr_av(k,j,i) = pr_av(k,j,i) + mean_r
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'pt' )
+             IF ( .NOT. cloud_physics ) THEN
+                DO  i = nxl-1, nxr+1
+                   DO  j = nys-1, nyn+1
+                      DO  k = nzb, nzt+1
+                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSE
+                DO  i = nxl-1, nxr+1
+                   DO  j = nys-1, nyn+1
+                      DO  k = nzb, nzt+1
+                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i) + l_d_cp * &
+                                                       pt_d_t(k) * ql(k,j,i)
+                      ENDDO
                    ENDDO
                 ENDDO
              ENDIF
 
-          CASE ( 'e' )
-             IF ( ALLOCATED( e_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         e_av(k,j,i) = e_av(k,j,i) + e(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'thetal' )
-             IF ( ALLOCATED( lpt_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         lpt_av(k,j,i) = lpt_av(k,j,i) + pt(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'lwp*' )
-             IF ( ALLOCATED( lwp_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i) * dzw(1:nzt+1) )            &
-                                                                       * rho_surface
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'ol*' )
-             IF ( ALLOCATED( ol_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         ol_av(j,i) = ol_av(j,i) + surf_def_h(0)%ol(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         ol_av(j,i) = ol_av(j,i) + surf_lsm_h(0)%ol(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         ol_av(j,i) = ol_av(j,i) + surf_usm_h(0)%ol(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'p' )
-             IF ( ALLOCATED( p_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         p_av(k,j,i) = p_av(k,j,i) + p(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'pc' )
-             IF ( ALLOCATED( pc_av ) ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
-                         pc_av(k,j,i) = pc_av(k,j,i) + prt_count(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'pr' )
-             IF ( ALLOCATED( pr_av ) ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
-                         number_of_particles = prt_count(k,j,i)
-                         IF ( number_of_particles <= 0 )  CYCLE
-                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
-                         s_r2 = 0.0_wp
-                         s_r3 = 0.0_wp
-
-                         DO  n = 1, number_of_particles
-                            IF ( particles(n)%particle_mask )  THEN
-                               s_r2 = s_r2 + particles(n)%radius**2 * particles(n)%weight_factor
-                               s_r3 = s_r3 + particles(n)%radius**3 * particles(n)%weight_factor
-                            ENDIF
-                         ENDDO
-
-                         IF ( s_r2 > 0.0_wp )  THEN
-                            mean_r = s_r3 / s_r2
-                         ELSE
-                            mean_r = 0.0_wp
-                         ENDIF
-                         pr_av(k,j,i) = pr_av(k,j,i) + mean_r
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'theta' )
-             IF ( ALLOCATED( pt_av ) ) THEN
-                IF ( .NOT. bulk_cloud_model ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                            pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
-                         ENDDO
-                      ENDDO
-                   ENDDO
-                ELSE
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                            pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i) + lv_d_cp * d_exner(k)         &
-                                                                              * ql(k,j,i)
-                         ENDDO
-                      ENDDO
-                   ENDDO
-                ENDIF
-             ENDIF
-
           CASE ( 'q' )
-             IF ( ALLOCATED( q_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+             
           CASE ( 'ql' )
-             IF ( ALLOCATED( ql_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
 
           CASE ( 'ql_c' )
-             IF ( ALLOCATED( ql_c_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         ql_c_av(k,j,i) = ql_c_av(k,j,i) + ql_c(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      ql_c_av(k,j,i) = ql_c_av(k,j,i) + ql_c(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
 
           CASE ( 'ql_v' )
-             IF ( ALLOCATED( ql_v_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         ql_v_av(k,j,i) = ql_v_av(k,j,i) + ql_v(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      ql_v_av(k,j,i) = ql_v_av(k,j,i) + ql_v(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
 
           CASE ( 'ql_vp' )
-             IF ( ALLOCATED( ql_vp_av ) ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
-                         number_of_particles = prt_count(k,j,i)
-                         IF ( number_of_particles <= 0 )  CYCLE
-                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
-                         DO  n = 1, number_of_particles
-                            IF ( particles(n)%particle_mask )  THEN
-                               ql_vp_av(k,j,i) = ql_vp_av(k,j,i) +                                 &
-                                                 particles(n)%weight_factor / number_of_particles
-                            ENDIF
-                         ENDDO
-                      ENDDO
-                   ENDDO
-                ENDDO
-             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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%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(0)%end_index(j,i)
-                         qsurf_av(j,i) = qsurf_av(j,i) + surf_lsm_h(0)%q_surface(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         qsurf_av(j,i) = qsurf_av(j,i) + surf_usm_h(0)%q_surface(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'qsws*' )
-!
-!--          In case of default surfaces, clean-up flux by density.
-!--          In case of land- and urban-surfaces, convert fluxes into dynamic units.
-!--          Question (maronga): are the .NOT. statements really required?
-             IF ( ALLOCATED( qsws_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         qsws_av(j,i) = qsws_av(j,i) + surf_def_h(0)%qsws(m) *                     &
-                                        waterflux_output_conversion(nzb)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         qsws_av(j,i) = qsws_av(j,i) + surf_lsm_h(0)%qsws(m) * l_v
-                      ELSEIF ( match_usm  .AND.  .NOT. match_lsm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         qsws_av(j,i) = qsws_av(j,i) + surf_usm_h(0)%qsws(m) * l_v
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + ql_vp(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
 
           CASE ( 'qv' )
-             IF ( ALLOCATED( qv_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         qv_av(k,j,i) = qv_av(k,j,i) + q(k,j,i) - ql(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'r_a*' )
-             IF ( ALLOCATED( r_a_av ) ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      match_lsm = surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         r_a_av(j,i) = r_a_av(j,i) + surf_lsm_h(0)%r_a(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         r_a_av(j,i) = r_a_av(j,i) + surf_usm_h(0)%frac(m,ind_veg_wall)  *            &
-                                                     surf_usm_h(0)%r_a(m)       +                     &
-                                                     surf_usm_h(0)%frac(m,ind_pav_green) *            &
-                                                     surf_usm_h(0)%r_a_green(m) +                     &
-                                                     surf_usm_h(0)%frac(m,ind_wat_win)   *            &
-                                                     surf_usm_h(0)%r_a_window(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      qv_av(k,j,i) = qv_av(k,j,i) + q(k,j,i) - ql(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
 
           CASE ( 's' )
-             IF ( ALLOCATED( s_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         s_av(k,j,i) = s_av(k,j,i) + s(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'shf*' )
-!
-!--          In case of default surfaces, clean-up flux by density.
-!--          In case of land- and urban-surfaces, convert fluxes into dynamic units.
-             IF ( ALLOCATED( shf_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         shf_av(j,i) = shf_av(j,i) + surf_def_h(0)%shf(m)  *                       &
-                                       heatflux_output_conversion(nzb)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         shf_av(j,i) = shf_av(j,i) + surf_lsm_h(0)%shf(m) * c_p
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         shf_av(j,i) = shf_av(j,i) + surf_usm_h(0)%shf(m) * c_p
-                      ENDIF
-                   ENDDO
-                ENDDO
-             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
-                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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         ssws_av(j,i) = ssws_av(j,i) + surf_def_h(0)%ssws(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         ssws_av(j,i) = ssws_av(j,i) + surf_lsm_h(0)%ssws(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         ssws_av(j,i) = ssws_av(j,i) + surf_usm_h(0)%ssws(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
-
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      s_av(k,j,i) = s_av(k,j,i) + q(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+             
           CASE ( 't*' )
-             IF ( ALLOCATED( ts_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         ts_av(j,i) = ts_av(j,i) + surf_def_h(0)%ts(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         ts_av(j,i) = ts_av(j,i) + surf_lsm_h(0)%ts(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         ts_av(j,i) = ts_av(j,i) + surf_usm_h(0)%ts(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'tsurf*' )
-             IF ( ALLOCATED( tsurf_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         tsurf_av(j,i) = tsurf_av(j,i) + surf_def_h(0)%pt_surface(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         tsurf_av(j,i) = tsurf_av(j,i) + surf_lsm_h(0)%pt_surface(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         tsurf_av(j,i) = tsurf_av(j,i) + surf_usm_h(0)%pt_surface(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   ts_av(j,i) = ts_av(j,i) + ts(j,i)
+                ENDDO
+             ENDDO
 
           CASE ( 'u' )
-             IF ( ALLOCATED( u_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'us*' )
-             IF ( ALLOCATED( us_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         us_av(j,i) = us_av(j,i) + surf_def_h(0)%us(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         us_av(j,i) = us_av(j,i) + surf_lsm_h(0)%us(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         us_av(j,i) = us_av(j,i) + surf_usm_h(0)%us(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'u*' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   us_av(j,i) = us_av(j,i) + us(j,i)
+                ENDDO
+             ENDDO
 
           CASE ( 'v' )
-             IF ( ALLOCATED( v_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'thetav' )
-             IF ( ALLOCATED( vpt_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         vpt_av(k,j,i) = vpt_av(k,j,i) + vpt(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'vpt' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      vpt_av(k,j,i) = vpt_av(k,j,i) + vpt(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
 
           CASE ( 'w' )
-             IF ( ALLOCATED( w_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'z0*' )
-             IF ( ALLOCATED( z0_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         z0_av(j,i) = z0_av(j,i) + surf_def_h(0)%z0(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         z0_av(j,i) = z0_av(j,i) + surf_lsm_h(0)%z0(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         z0_av(j,i) = z0_av(j,i) + surf_usm_h(0)%z0(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'z0h*' )
-             IF ( ALLOCATED( z0h_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         z0h_av(j,i) = z0h_av(j,i) + surf_def_h(0)%z0h(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         z0h_av(j,i) = z0h_av(j,i) + surf_lsm_h(0)%z0h(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         z0h_av(j,i) = z0h_av(j,i) + surf_usm_h(0)%z0h(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'z0q*' )
-             IF ( ALLOCATED( z0q_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(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
-                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         z0q_av(j,i) = z0q_av(j,i) + surf_def_h(0)%z0q(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h(0)%end_index(j,i)
-                         z0q_av(j,i) = z0q_av(j,i) + surf_lsm_h(0)%z0q(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h(0)%end_index(j,i)
-                         z0q_av(j,i) = z0q_av(j,i) + surf_usm_h(0)%z0q(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
 
           CASE DEFAULT
-
-!--          In case of urban surface variables it should be always checked if respective arrays are
-!--          allocated, at least in case of a restart run, as averaged usm arrays are not read from
-!--          file at the moment.
-             IF ( urban_surface )  THEN
-                CALL usm_3d_data_averaging( 'allocate', trimvar )
-             ENDIF
-
-!
-!--          Summing up data from all other modules
-             CALL module_interface_3d_data_averaging( 'sum', trimvar )
-
+!
+!--          User-defined quantity
+             CALL user_3d_data_averaging( 'sum', doav(ii) )
 
        END SELECT
@@ -1003,5 +379,5 @@
     ENDDO
 
-    CALL cpu_log( log_point(34), 'sum_up_3d_data', 'stop' )
+    CALL cpu_log (log_point(34),'sum_up_3d_data','stop','nobarrier')
 
 
