Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4516)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4591)
@@ -1,64 +1,66 @@
 !> @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/>.
+! 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-2020 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+!
 !
 ! Current revisions:
-! ------------------
-! 
-! 
+! -----------------
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
-! remove double index
-! 
+! 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.
-! 
+! 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
-! 
+! 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
@@ -70,52 +72,114 @@
 ! Description:
 ! ------------
-!> Sum-up the values of 3d-arrays. The real averaging is later done in routine
-!> average_3d_data. 
-!------------------------------------------------------------------------------!
+!> 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,                         &
+
+
+    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,                                                  &
+    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,     &
+    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,           &
+    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,                                                      &
+    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,                                                     &
+    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
 
@@ -123,8 +187,6 @@
     IMPLICIT NONE
 
-    LOGICAL      ::  match_def !< flag indicating default-type surface
-    LOGICAL      ::  match_lsm !< flag indicating natural-type surface
-    LOGICAL      ::  match_usm !< flag indicating urban-type surface
-    
+    CHARACTER(LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
+
     INTEGER(iwp) ::  i   !< grid index x direction
     INTEGER(iwp) ::  ii  !< running index
@@ -134,9 +196,12 @@
     INTEGER(iwp) ::  n   !< running index over number of particles per grid box
 
-    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
-
-    CHARACTER (LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
+    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
+
 
 
@@ -144,7 +209,6 @@
 
 !
-!-- 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
+!-- 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)
     IF ( average_count_3d == 0 )  THEN
@@ -286,5 +350,5 @@
                    ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
                 ENDIF
-                ssws_av = 0.0_wp                
+                ssws_av = 0.0_wp
 
              CASE ( 't*' )
@@ -376,29 +440,23 @@
                    DO  j = nys, nyn
 !
-!--                   Check whether grid point is a natural- or urban-type 
-!--                   surface.
-                      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)
-!
-!--                   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
+!--                   Check whether grid point is a natural- or urban-type surface.
+                      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)
+!
+!--                   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%end_index(j,i)
-                         ghf_av(j,i) = ghf_av(j,i) +                           &
-                                         surf_lsm_h%ghf(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         ghf_av(j,i) = ghf_av(j,i) +                           &
-                                         surf_usm_h%frac(m,ind_veg_wall)  *    &
-                                         surf_usm_h%wghf_eb(m)        +        &
-                                         surf_usm_h%frac(m,ind_pav_green) *    &
-                                         surf_usm_h%wghf_eb_green(m)  +        &
-                                         surf_usm_h%frac(m,ind_wat_win)   *    &
-                                         surf_usm_h%wghf_eb_window(m)
+                         ghf_av(j,i) = ghf_av(j,i) + surf_lsm_h%ghf(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         ghf_av(j,i) = ghf_av(j,i) + surf_usm_h%frac(m,ind_veg_wall)  *            &
+                                                     surf_usm_h%wghf_eb(m)        +                &
+                                                     surf_usm_h%frac(m,ind_pav_green) *            &
+                                                     surf_usm_h%wghf_eb_green(m)  +                &
+                                                     surf_usm_h%frac(m,ind_wat_win)   *            &
+                                                     surf_usm_h%wghf_eb_window(m)
                       ENDIF
                    ENDDO
@@ -432,6 +490,6 @@
                 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
+                      lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i) * dzw(1:nzt+1) )            &
+                                                                       * rho_surface
                    ENDDO
                 ENDDO
@@ -442,23 +500,17 @@
                 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)
+                      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)
-                         ol_av(j,i) = ol_av(j,i) +                             &
-                                         surf_def_h(0)%ol(m)
+                         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%end_index(j,i)
-                         ol_av(j,i) = ol_av(j,i) +                             &
-                                         surf_lsm_h%ol(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         ol_av(j,i) = ol_av(j,i) +                             &
-                                         surf_usm_h%ol(m)
+                         ol_av(j,i) = ol_av(j,i) + surf_lsm_h%ol(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         ol_av(j,i) = ol_av(j,i) + surf_usm_h%ol(m)
                       ENDIF
                    ENDDO
@@ -495,6 +547,5 @@
                          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)
+                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
                          s_r2 = 0.0_wp
                          s_r3 = 0.0_wp
@@ -502,8 +553,6 @@
                          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
+                               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
@@ -534,6 +583,6 @@
                    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)
+                            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
@@ -565,5 +614,5 @@
 
           CASE ( 'ql_c' )
-             IF ( ALLOCATED( ql_c_av ) ) THEN 
+             IF ( ALLOCATED( ql_c_av ) ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -576,5 +625,5 @@
 
           CASE ( 'ql_v' )
-             IF ( ALLOCATED( ql_v_av ) ) THEN 
+             IF ( ALLOCATED( ql_v_av ) ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -587,5 +636,5 @@
 
           CASE ( 'ql_vp' )
-             IF ( ALLOCATED( ql_vp_av ) ) THEN 
+             IF ( ALLOCATED( ql_vp_av ) ) THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -593,11 +642,9 @@
                          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)
+                         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
+                               ql_vp_av(k,j,i) = ql_vp_av(k,j,i) +                                 &
+                                                 particles(n)%weight_factor / number_of_particles
                             ENDIF
                          ENDDO
@@ -611,23 +658,17 @@
                 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)
+                      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)
+                         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)
+                         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
@@ -638,30 +679,23 @@
 !
 !--          In case of default surfaces, clean-up flux by density.
-!--          In case of land- and urban-surfaces, convert fluxes into
-!--          dynamic units.
+!--          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%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 ( 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%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)
-                         qsws_av(j,i) = qsws_av(j,i) +                         &
-                                         surf_def_h(0)%qsws(m) *               &
-                                         waterflux_output_conversion(nzb)
+                         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%end_index(j,i)
-                         qsws_av(j,i) = qsws_av(j,i) +                         &
-                                         surf_lsm_h%qsws(m) * l_v
+                         qsws_av(j,i) = qsws_av(j,i) + surf_lsm_h%qsws(m) * l_v
                       ELSEIF ( match_usm  .AND.  .NOT. match_lsm )  THEN
                          m = surf_usm_h%end_index(j,i)
-                         qsws_av(j,i) = qsws_av(j,i) +                         &
-                                         surf_usm_h%qsws(m) * l_v
+                         qsws_av(j,i) = qsws_av(j,i) + surf_usm_h%qsws(m) * l_v
                       ENDIF
                    ENDDO
@@ -670,5 +704,5 @@
 
           CASE ( 'qv' )
-             IF ( ALLOCATED( qv_av ) ) THEN 
+             IF ( ALLOCATED( qv_av ) ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -681,25 +715,21 @@
 
           CASE ( 'r_a*' )
-             IF ( ALLOCATED( r_a_av ) ) THEN 
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      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 ( ALLOCATED( r_a_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      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_lsm  .AND.  .NOT. match_usm )  THEN
                          m = surf_lsm_h%end_index(j,i)
-                         r_a_av(j,i) = r_a_av(j,i) +                           &
-                                         surf_lsm_h%r_a(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         r_a_av(j,i) = r_a_av(j,i) +                           &
-                                         surf_usm_h%frac(m,ind_veg_wall)  *    &
-                                         surf_usm_h%r_a(m)       +             & 
-                                         surf_usm_h%frac(m,ind_pav_green) *    &
-                                         surf_usm_h%r_a_green(m) +             & 
-                                         surf_usm_h%frac(m,ind_wat_win)   *    &
-                                         surf_usm_h%r_a_window(m)
+                         r_a_av(j,i) = r_a_av(j,i) + surf_lsm_h%r_a(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         r_a_av(j,i) = r_a_av(j,i) + surf_usm_h%frac(m,ind_veg_wall)  *            &
+                                                     surf_usm_h%r_a(m)       +                     &
+                                                     surf_usm_h%frac(m,ind_pav_green) *            &
+                                                     surf_usm_h%r_a_green(m) +                     &
+                                                     surf_usm_h%frac(m,ind_wat_win)   *            &
+                                                     surf_usm_h%r_a_window(m)
                       ENDIF
                    ENDDO
@@ -708,5 +738,5 @@
 
           CASE ( 's' )
-             IF ( ALLOCATED( s_av ) ) THEN 
+             IF ( ALLOCATED( s_av ) ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -721,29 +751,22 @@
 !
 !--          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%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)
+!--          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%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)
-                         shf_av(j,i) = shf_av(j,i) +                           &
-                                         surf_def_h(0)%shf(m)  *               &
-                                         heatflux_output_conversion(nzb)
+                         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%end_index(j,i)
-                         shf_av(j,i) = shf_av(j,i) +                           &
-                                         surf_lsm_h%shf(m) * c_p
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         shf_av(j,i) = shf_av(j,i) +                           &
-                                         surf_usm_h%shf(m) * c_p
+                         shf_av(j,i) = shf_av(j,i) + surf_lsm_h%shf(m) * c_p
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         shf_av(j,i) = shf_av(j,i) + surf_usm_h%shf(m) * c_p
                       ENDIF
                    ENDDO
@@ -762,26 +785,20 @@
 
           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%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 ( 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%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)
-                         ssws_av(j,i) = ssws_av(j,i) +                         &
-                                         surf_def_h(0)%ssws(m)
+                         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%end_index(j,i)
-                         ssws_av(j,i) = ssws_av(j,i) +                         &
-                                         surf_lsm_h%ssws(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         ssws_av(j,i) = ssws_av(j,i) +                         &
-                                         surf_usm_h%ssws(m)
+                         ssws_av(j,i) = ssws_av(j,i) + surf_lsm_h%ssws(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         ssws_av(j,i) = ssws_av(j,i) + surf_usm_h%ssws(m)
                       ENDIF
                    ENDDO
@@ -790,26 +807,20 @@
 
           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%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 ( 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%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)
-                         ts_av(j,i) = ts_av(j,i) +                             &
-                                         surf_def_h(0)%ts(m)
+                         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%end_index(j,i)
-                         ts_av(j,i) = ts_av(j,i) +                             &
-                                         surf_lsm_h%ts(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         ts_av(j,i) = ts_av(j,i) +                             &
-                                         surf_usm_h%ts(m)
+                         ts_av(j,i) = ts_av(j,i) + surf_lsm_h%ts(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         ts_av(j,i) = ts_av(j,i) + surf_usm_h%ts(m)
                       ENDIF
                    ENDDO
@@ -818,26 +829,20 @@
 
           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%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 ( 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%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)
-                         tsurf_av(j,i) = tsurf_av(j,i) +                       &
-                                         surf_def_h(0)%pt_surface(m)
+                         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%end_index(j,i)
-                         tsurf_av(j,i) = tsurf_av(j,i) +                       &
-                                         surf_lsm_h%pt_surface(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         tsurf_av(j,i) = tsurf_av(j,i) +                       &
-                                         surf_usm_h%pt_surface(m)
+                         tsurf_av(j,i) = tsurf_av(j,i) + surf_lsm_h%pt_surface(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         tsurf_av(j,i) = tsurf_av(j,i) + surf_usm_h%pt_surface(m)
                       ENDIF
                    ENDDO
@@ -857,26 +862,20 @@
 
           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%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 ( 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%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)
-                         us_av(j,i) = us_av(j,i) +                             &
-                                         surf_def_h(0)%us(m)
+                         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%end_index(j,i)
-                         us_av(j,i) = us_av(j,i) +                             &
-                                         surf_lsm_h%us(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         us_av(j,i) = us_av(j,i) +                             &
-                                         surf_usm_h%us(m)
+                         us_av(j,i) = us_av(j,i) + surf_lsm_h%us(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         us_av(j,i) = us_av(j,i) + surf_usm_h%us(m)
                       ENDIF
                    ENDDO
@@ -885,5 +884,5 @@
 
           CASE ( 'v' )
-             IF ( ALLOCATED( v_av ) ) THEN 
+             IF ( ALLOCATED( v_av ) ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -896,5 +895,5 @@
 
           CASE ( 'thetav' )
-             IF ( ALLOCATED( vpt_av ) ) THEN 
+             IF ( ALLOCATED( vpt_av ) ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -907,5 +906,5 @@
 
           CASE ( 'w' )
-             IF ( ALLOCATED( w_av ) ) THEN 
+             IF ( ALLOCATED( w_av ) ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -918,82 +917,64 @@
 
           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%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 ( 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%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)
-                         z0_av(j,i) = z0_av(j,i) +                             &
-                                         surf_def_h(0)%z0(m)
+                         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%end_index(j,i)
-                         z0_av(j,i) = z0_av(j,i) +                             &
-                                         surf_lsm_h%z0(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         z0_av(j,i) = z0_av(j,i) +                             &
-                                         surf_usm_h%z0(m)
-                      ENDIF
-                   ENDDO
-                ENDDO   
+                         z0_av(j,i) = z0_av(j,i) + surf_lsm_h%z0(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         z0_av(j,i) = z0_av(j,i) + surf_usm_h%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%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 ( 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%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)
-                         z0h_av(j,i) = z0h_av(j,i) +                           &
-                                         surf_def_h(0)%z0h(m)
+                         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%end_index(j,i)
-                         z0h_av(j,i) = z0h_av(j,i) +                           &
-                                         surf_lsm_h%z0h(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         z0h_av(j,i) = z0h_av(j,i) +                           &
-                                         surf_usm_h%z0h(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
-   
+                         z0h_av(j,i) = z0h_av(j,i) + surf_lsm_h%z0h(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         z0h_av(j,i) = z0h_av(j,i) + surf_usm_h%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%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 ( 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%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)
-                         z0q_av(j,i) = z0q_av(j,i) +                           &
-                                         surf_def_h(0)%z0q(m)
+                         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%end_index(j,i)
-                         z0q_av(j,i) = z0q_av(j,i) +                           &
-                                         surf_lsm_h%z0q(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         z0q_av(j,i) = z0q_av(j,i) +                           &
-                                         surf_usm_h%z0q(m)
+                         z0q_av(j,i) = z0q_av(j,i) + surf_lsm_h%z0q(m)
+                      ELSEIF ( match_usm )  THEN
+                         m = surf_usm_h%end_index(j,i)
+                         z0q_av(j,i) = z0q_av(j,i) + surf_usm_h%z0q(m)
                       ENDIF
                    ENDDO
@@ -1003,7 +984,7 @@
           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.
+!--          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 )
