Index: palm/trunk/SOURCE/data_output_tseries.f90
===================================================================
--- palm/trunk/SOURCE/data_output_tseries.f90	(revision 4582)
+++ palm/trunk/SOURCE/data_output_tseries.f90	(revision 4583)
@@ -1,20 +1,19 @@
 !> @file data_output_tseries.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 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.
+! 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/>.
+! 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:
@@ -25,6 +24,9 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! unused format removed
@@ -36,15 +38,15 @@
 ! Description:
 ! ------------
-!> Time series output for PROFIL. Always all time series are stored. A selection
-!> can be applied via the PROFIL-parameters in close_file.
-!------------------------------------------------------------------------------!
+!> Time series output for PROFIL. Always all time series are stored. A selection can be applied via
+!> the PROFIL-parameters in close_file.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE data_output_tseries
- 
 
-    USE control_parameters,                                                    &
+
+    USE control_parameters,                                                                        &
         ONLY:  dots_time_count, time_since_reference_point
 
-    USE cpulog,                                                                &
-        ONLY:  cpu_log, log_point  
+    USE cpulog,                                                                                    &
+        ONLY:  cpu_log, log_point
 
     USE kinds
@@ -53,13 +55,12 @@
     USE NETCDF
 #endif
-    USE netcdf_interface,                                                      &
-        ONLY:  dots_num, id_set_ts, id_var_dots, id_var_time_ts, nc_stat,      &
-               netcdf_handle_error
+    USE netcdf_interface,                                                                          &
+        ONLY:  dots_num, id_set_ts, id_var_dots, id_var_time_ts, nc_stat, netcdf_handle_error
 
     USE pegrid
 
     USE profil_parameter
-    
-    USE statistics,                                                            &
+
+    USE statistics,                                                                                &
         ONLY:  flow_statistics_called, statistic_regions, ts_value
 
@@ -84,17 +85,16 @@
 !--    Open file for time series output in NetCDF format
        CALL check_open( 105 )
-       
+
 !--    Increment the counter for number of output times
-!      CAUTION: The following line has to be after the call of the subroutine
-!               check_open, since check_open resets the counter dots_time_count
-!               to 0, if a new file is opened
+!--    CAUTION: The following line has to be after the call of the subroutine check_open, since
+!--             check_open resets the counter dots_time_count to 0, if a new file is opened
        dots_time_count = dots_time_count + 1
-       
+
 #if defined( __netcdf )
 !
 !--    Update the time series time axis
-       nc_stat = NF90_PUT_VAR( id_set_ts, id_var_time_ts,        &
-                               (/ time_since_reference_point /), &
-                               start = (/ dots_time_count /),    &
+       nc_stat = NF90_PUT_VAR( id_set_ts, id_var_time_ts,                                          &
+                               (/ time_since_reference_point /),                                   &
+                               start = (/ dots_time_count /),                                      &
                                count = (/ 1 /) )
        CALL netcdf_handle_error( 'data_output_tseries', 350 )
@@ -108,7 +108,7 @@
 #if defined( __netcdf )
           DO  i = 1, dots_num
-             nc_stat = NF90_PUT_VAR( id_set_ts, id_var_dots(i,sr),  &
-                                     (/ ts_value(i,sr) /),          &
-                                     start = (/ dots_time_count /), &
+             nc_stat = NF90_PUT_VAR( id_set_ts, id_var_dots(i,sr),                                 &
+                                     (/ ts_value(i,sr) /),                                         &
+                                     start = (/ dots_time_count /),                                &
                                      count = (/ 1 /) )
              CALL netcdf_handle_error( 'data_output_tseries', 351 )
Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4582)
+++ palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4583)
@@ -1,20 +1,19 @@
 !> @file diagnostic_output_quantities_mod.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:
@@ -25,30 +24,32 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4560 2020-06-11 12:19:47Z suehring
 ! - Bugfix in calculation of vertical momentum and scalar fluxes
 ! - remove averaged output variables from PUBLIC list
-! 
+!
 ! 4535 2020-05-15 12:07:23Z raasch
 ! bugfix for restart data format query
-! 
+!
 ! 4518 2020-05-04 15:44:28Z suehring
-! * Define arrays over ghost points in order to allow for standard mpi-io 
-!   treatment. By this modularization of restart-data input is possible with
-!   the module interface.
+! * Define arrays over ghost points in order to allow for standard mpi-io treatment. By this
+!   modularization of restart-data input is possible with the module interface.
 ! * Bugfix: add missing restart input of wtheta_av, wq_av, wu_av, and wv_av.
-! 
+!
 ! 4517 2020-05-03 14:29:30Z raasch
 ! use statement for exchange horiz added,
 ! bugfix for call of exchange horiz 2d
-! 
+!
 ! 4431 2020-02-27 23:23:01Z gronemeier
 ! added wspeed and wdir output; bugfix: set fill_value in case of masked output
 !
 ! 4360 2020-01-07 11:25:50Z suehring
-! added output of wu, wv, wtheta and wq to enable covariance calculation
-! according to temporal EC method
+! added output of wu, wv, wtheta and wq to enable covariance calculation according to temporal EC
+! method
 !
 ! 4346 2019-12-18 11:55:56Z motisi
-! Introduction of wall_flags_total_0, which currently sets bits based on static
-! topography information used in wall_flags_static_0
+! Introduction of wall_flags_total_0, which currently sets bits based on static topography
+! information used in wall_flags_static_0
 !
 ! 4331 2019-12-10 18:25:02Z suehring
@@ -63,10 +64,9 @@
 !
 ! 4167 2019-08-16 11:01:48Z suehring
-! Changed behaviour of masked output over surface to follow terrain and ignore
-! buildings (J.Resler, T.Gronemeier)
+! Changed behaviour of masked output over surface to follow terrain and ignore buildings
+! (J.Resler, T.Gronemeier)
 !
 ! 4157 2019-08-14 09:19:12Z suehring
-! Initialization restructured, in order to work also when data output during
-! spin-up is enabled.
+! Initialization restructured, in order to work also when data output during spin-up is enabled.
 !
 ! 4132 2019-08-02 12:34:17Z suehring
@@ -74,10 +74,9 @@
 !
 ! 4069 2019-07-01 14:05:51Z Giersch
-! Masked output running index mid has been introduced as a local variable to
-! avoid runtime error (Loop variable has been modified) in time_integration
+! Masked output running index mid has been introduced as a local variable to avoid runtime error
+! (Loop variable has been modified) in time_integration
 !
 ! 4039 2019-06-18 10:32:41Z suehring
-! - Add output of uu, vv, ww to enable variance calculation according temporal
-!   EC method
+! - Add output of uu, vv, ww to enable variance calculation according temporal EC method
 ! - Allocate arrays only when they are required
 ! - Formatting adjustment
@@ -89,6 +88,6 @@
 !
 ! 3995 2019-05-22 18:59:54Z suehring
-! Avoid compiler warnings about unused variable and fix string operation which
-! is not allowed with PGI compiler
+! Avoid compiler warnings about unused variable and fix string operation which is not allowed with
+! PGI compiler
 !
 ! 3994 2019-05-22 18:08:09Z suehring
@@ -103,48 +102,48 @@
 ! ------------
 !> ...
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  MODULE diagnostic_output_quantities_mod
 
-    USE arrays_3d,                                                             &
-        ONLY:  ddzu,                                                           &
-               pt,                                                             &
-               q,                                                              &
-               u,                                                              &
-               v,                                                              &
-               w,                                                              &
-               zu,                                                             &
+    USE arrays_3d,                                                                                 &
+        ONLY:  ddzu,                                                                               &
+               pt,                                                                                 &
+               q,                                                                                  &
+               u,                                                                                  &
+               v,                                                                                  &
+               w,                                                                                  &
+               zu,                                                                                 &
                zw
 
-    USE basic_constants_and_equations_mod,                                     &
+    USE basic_constants_and_equations_mod,                                                         &
         ONLY:  kappa, pi
 
-    USE control_parameters,                                                    &
-        ONLY:  current_timestep_number,                                        &
-               data_output,                                                    &
-               message_string,                                                 &
-               restart_data_format_output,                                     &
+    USE control_parameters,                                                                        &
+        ONLY:  current_timestep_number,                                                            &
+               data_output,                                                                        &
+               message_string,                                                                     &
+               restart_data_format_output,                                                         &
                varnamelength
 !
-!     USE cpulog,                                                                &
+!     USE cpulog,                                                                                   &
 !         ONLY:  cpu_log, log_point
 
-    USE exchange_horiz_mod,                                                    &
+    USE exchange_horiz_mod,                                                                        &
         ONLY:  exchange_horiz_2d
 
-   USE grid_variables,                                                         &
+   USE grid_variables,                                                                             &
         ONLY:  ddx, ddy
 
-    USE indices,                                                               &
-        ONLY:  nbgp,                                                           &
-               nxl,                                                            &
-               nxlg,                                                           &
-               nxr,                                                            &
-               nxrg,                                                           &
-               nyn,                                                            &
-               nyng,                                                           &
-               nys,                                                            &
-               nysg,                                                           &
-               nzb,                                                            &
-               nzt,                                                            &
+    USE indices,                                                                                   &
+        ONLY:  nbgp,                                                                               &
+               nxl,                                                                                &
+               nxlg,                                                                               &
+               nxr,                                                                                &
+               nxrg,                                                                               &
+               nyn,                                                                                &
+               nyng,                                                                               &
+               nys,                                                                                &
+               nysg,                                                                               &
+               nzb,                                                                                &
+               nzt,                                                                                &
                wall_flags_total_0
 
@@ -156,8 +155,8 @@
                wrd_mpi_io
 
-    USE surface_mod,                                                           &
-        ONLY:  surf_def_h,                                                     &
-               surf_lsm_h,                                                     &
-               surf_type,                                                      &
+    USE surface_mod,                                                                               &
+        ONLY:  surf_def_h,                                                                         &
+               surf_lsm_h,                                                                         &
+               surf_type,                                                                          &
                surf_usm_h
 
@@ -209,19 +208,19 @@
 !
 !-- Public variables
-    PUBLIC do_all,                                                             &
-           initialized_diagnostic_output_quantities,                           &
-           prepared_diagnostic_output_quantities,                              &
+    PUBLIC do_all,                                                                                 &
+           initialized_diagnostic_output_quantities,                                               &
+           prepared_diagnostic_output_quantities,                                                  &
            timestep_number_at_prev_calc
 !
 !-- Public routines
-    PUBLIC doq_3d_data_averaging,                                              &
-           doq_calculate,                                                      &
-           doq_check_data_output,                                              &
-           doq_define_netcdf_grid,                                             &
-           doq_init,                                                           &
-           doq_output_2d,                                                      &
-           doq_output_3d,                                                      &
-           doq_output_mask,                                                    &
-           doq_rrd_local,                                                      &
+    PUBLIC doq_3d_data_averaging,                                                                  &
+           doq_calculate,                                                                          &
+           doq_check_data_output,                                                                  &
+           doq_define_netcdf_grid,                                                                 &
+           doq_init,                                                                               &
+           doq_output_2d,                                                                          &
+           doq_output_3d,                                                                          &
+           doq_output_mask,                                                                        &
+           doq_rrd_local,                                                                          &
            doq_wrd_local
 
@@ -275,13 +274,13 @@
  CONTAINS
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Sum up and time-average diagnostic output quantities as well as allocate
-!> the array necessary for storing the average.
-!------------------------------------------------------------------------------!
+!> Sum up and time-average diagnostic output quantities as well as allocate the array necessary for
+!> storing the average.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE doq_3d_data_averaging( mode, variable )
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  average_count_3d
 
@@ -293,4 +292,5 @@
     INTEGER(iwp) ::  k !<
 
+
     IF ( mode == 'allocate' )  THEN
 
@@ -383,5 +383,5 @@
 
           CASE ( 'ti' )
-             IF ( ALLOCATED( ti_av ) ) THEN
+             IF ( ALLOCATED( ti_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -394,5 +394,5 @@
 
           CASE ( 'uu' )
-             IF ( ALLOCATED( uu_av ) ) THEN
+             IF ( ALLOCATED( uu_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -405,5 +405,5 @@
 
           CASE ( 'vv' )
-             IF ( ALLOCATED( vv_av ) ) THEN
+             IF ( ALLOCATED( vv_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -416,5 +416,5 @@
 
           CASE ( 'ww' )
-             IF ( ALLOCATED( ww_av ) ) THEN
+             IF ( ALLOCATED( ww_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -427,5 +427,5 @@
 
           CASE ( 'wu' )
-             IF ( ALLOCATED( wu_av ) ) THEN
+             IF ( ALLOCATED( wu_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -438,5 +438,5 @@
 
           CASE ( 'wv' )
-             IF ( ALLOCATED( wv_av ) ) THEN
+             IF ( ALLOCATED( wv_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -449,5 +449,5 @@
 
           CASE ( 'wtheta' )
-             IF ( ALLOCATED( wtheta_av ) ) THEN
+             IF ( ALLOCATED( wtheta_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -460,5 +460,5 @@
 
           CASE ( 'wq' )
-             IF ( ALLOCATED( wq_av ) ) THEN
+             IF ( ALLOCATED( wq_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -471,5 +471,5 @@
 
           CASE ( 'theta_2m*' )
-             IF ( ALLOCATED( pt_2m_av ) ) THEN
+             IF ( ALLOCATED( pt_2m_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -480,5 +480,5 @@
 
           CASE ( 'wspeed_10m*' )
-             IF ( ALLOCATED( uv_10m_av ) ) THEN
+             IF ( ALLOCATED( uv_10m_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -489,5 +489,5 @@
 
           CASE ( 'wspeed' )
-            IF ( ALLOCATED( wspeed_av ) ) THEN
+            IF ( ALLOCATED( wspeed_av ) )  THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
@@ -500,5 +500,5 @@
 
           CASE ( 'wdir' )
-             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) ) THEN
+             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -521,5 +521,5 @@
 
           CASE ( 'ti' )
-             IF ( ALLOCATED( ti_av ) ) THEN
+             IF ( ALLOCATED( ti_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -532,5 +532,5 @@
 
           CASE ( 'uu' )
-             IF ( ALLOCATED( uu_av ) ) THEN
+             IF ( ALLOCATED( uu_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -543,5 +543,5 @@
 
           CASE ( 'vv' )
-             IF ( ALLOCATED( vv_av ) ) THEN
+             IF ( ALLOCATED( vv_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -554,5 +554,5 @@
 
           CASE ( 'ww' )
-             IF ( ALLOCATED( ww_av ) ) THEN
+             IF ( ALLOCATED( ww_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -565,5 +565,5 @@
 
           CASE ( 'wu' )
-             IF ( ALLOCATED( wu_av ) ) THEN
+             IF ( ALLOCATED( wu_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -576,5 +576,5 @@
 
           CASE ( 'wv' )
-             IF ( ALLOCATED( wv_av ) ) THEN
+             IF ( ALLOCATED( wv_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -587,5 +587,5 @@
 
           CASE ( 'wtheta' )
-             IF ( ALLOCATED( wtheta_av ) ) THEN
+             IF ( ALLOCATED( wtheta_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -598,5 +598,5 @@
 
           CASE ( 'wq' )
-             IF ( ALLOCATED( wq_av ) ) THEN
+             IF ( ALLOCATED( wq_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -609,5 +609,5 @@
 
          CASE ( 'theta_2m*' )
-            IF ( ALLOCATED( pt_2m_av ) ) THEN
+            IF ( ALLOCATED( pt_2m_av ) )  THEN
                DO  i = nxlg, nxrg
                   DO  j = nysg, nyng
@@ -619,5 +619,5 @@
 
          CASE ( 'wspeed_10m*' )
-            IF ( ALLOCATED( uv_10m_av ) ) THEN
+            IF ( ALLOCATED( uv_10m_av ) )  THEN
                DO  i = nxlg, nxrg
                   DO  j = nysg, nyng
@@ -629,5 +629,5 @@
 
          CASE ( 'wspeed' )
-             IF ( ALLOCATED( wspeed_av ) ) THEN
+             IF ( ALLOCATED( wspeed_av ) )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -640,5 +640,5 @@
 
           CASE ( 'wdir' )
-             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) ) THEN
+             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) )  THEN
 
                 IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
@@ -652,6 +652,6 @@
                          u_center_av(k,j,i) = u_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                          v_center_av(k,j,i) = v_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                         wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) ) &
-                                        / pi * 180.0_wp + 180.0_wp
+                         wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) )          &
+                                          / pi * 180.0_wp + 180.0_wp
                       ENDDO
                    ENDDO
@@ -666,9 +666,9 @@
  END SUBROUTINE doq_3d_data_averaging
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check data output for diagnostic output
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
 
@@ -682,4 +682,5 @@
     INTEGER(iwp), OPTIONAL, INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
 
+
     SELECT CASE ( TRIM( var ) )
 
@@ -719,6 +720,6 @@
 !--       Check if output quantity is _xy only.
           IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
-             message_string = 'illegal value for data_output: "' //            &
-                              TRIM( var ) // '" & only 2d-horizontal ' //      &
+             message_string = 'illegal value for data_output: "' //                                &
+                              TRIM( var ) // '" & only 2d-horizontal ' //                          &
                               'cross sections are allowed for this value'
              CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
@@ -736,10 +737,10 @@
  END SUBROUTINE doq_check_data_output
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !
 ! Description:
 ! ------------
 !> Subroutine defining appropriate grid for netcdf variables.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
 
@@ -747,9 +748,11 @@
 
     CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
-    LOGICAL, INTENT(OUT)           ::  found       !<
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
 
+    LOGICAL, INTENT(OUT)           ::  found       !<
+
+
     found  = .TRUE.
 
@@ -757,6 +760,6 @@
 !
 !--    s grid
-       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                 &
-              'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                 &
+       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                                     &
+              'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                                     &
               'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' )
 
@@ -787,9 +790,9 @@
 !
 !--    w grid
-       CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz',                                 &
-              'wu', 'wu_xy', 'wu_xz', 'wu_yz',                                 &
-              'wv', 'wv_xy', 'wv_xz', 'wv_yz',                                 &
-              'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz',                 &
-              'wq', 'wq_xy', 'wq_xz', 'wq_yz'  )
+       CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz',                                                     &
+              'wu', 'wu_xy', 'wu_xz', 'wu_yz',                                                     &
+              'wv', 'wv_xy', 'wv_xz', 'wv_yz',                                                     &
+              'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz',                                     &
+              'wq', 'wq_xy', 'wq_xz', 'wq_yz' )
 
           grid_x = 'x'
@@ -808,12 +811,12 @@
  END SUBROUTINE doq_define_netcdf_grid
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !
 ! Description:
 ! ------------
 !> Subroutine defining 2D output variables
-!------------------------------------------------------------------------------!
- SUBROUTINE doq_output_2d( av, variable, found, grid,                          &
-                           mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE doq_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do,       &
+                           fill_value )
 
 
@@ -841,4 +844,5 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
 
+
     flag_nr  = 0
     found    = .TRUE.
@@ -852,5 +856,5 @@
               to_be_resorted => ti
            ELSE
-              IF ( .NOT. ALLOCATED( ti_av ) ) THEN
+              IF ( .NOT. ALLOCATED( ti_av ) )  THEN
                  ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                  ti_av = REAL( fill_value, KIND = wp )
@@ -866,5 +870,5 @@
              to_be_resorted => uu
           ELSE
-             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
+             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
                 ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 uu_av = REAL( fill_value, KIND = wp )
@@ -880,5 +884,5 @@
              to_be_resorted => vv
           ELSE
-             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
+             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
                 ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 vv_av = REAL( fill_value, KIND = wp )
@@ -894,5 +898,5 @@
              to_be_resorted => ww
           ELSE
-             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
+             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
                 ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 ww_av = REAL( fill_value, KIND = wp )
@@ -908,5 +912,5 @@
              to_be_resorted => wu
           ELSE
-             IF ( .NOT. ALLOCATED( wu_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wu_av ) )  THEN
                 ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wu_av = REAL( fill_value, KIND = wp )
@@ -922,5 +926,5 @@
              to_be_resorted => wv
           ELSE
-             IF ( .NOT. ALLOCATED( wv_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wv_av ) )  THEN
                 ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wv_av = REAL( fill_value, KIND = wp )
@@ -936,5 +940,5 @@
              to_be_resorted => wtheta
           ELSE
-             IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
                 ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wtheta_av = REAL( fill_value, KIND = wp )
@@ -950,5 +954,5 @@
              to_be_resorted => wq
           ELSE
-             IF ( .NOT. ALLOCATED( wq_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wq_av ) )  THEN
                 ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wq_av = REAL( fill_value, KIND = wp )
@@ -968,5 +972,5 @@
              ENDDO
           ELSE
-             IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
+             IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
                 ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
                 pt_2m_av = REAL( fill_value, KIND = wp )
@@ -990,5 +994,5 @@
              ENDDO
           ELSE
-             IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
+             IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
                 ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
                 uv_10m_av = REAL( fill_value, KIND = wp )
@@ -1008,5 +1012,5 @@
              to_be_resorted => wspeed
           ELSE
-             IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
                 ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wspeed_av = REAL( fill_value, KIND = wp )
@@ -1022,5 +1026,5 @@
              to_be_resorted => wdir
           ELSE
-             IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
                 ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wdir_av = REAL( fill_value, KIND = wp )
@@ -1042,7 +1046,7 @@
           DO  j = nys, nyn
              DO  k = nzb_do, nzt_do
-                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
-                                     REAL( fill_value, KIND = wp ),            &
-                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
+                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                                    &
+                                         REAL( fill_value, KIND = wp ),                            &
+                                         BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
              ENDDO
           ENDDO
@@ -1053,12 +1057,11 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !
 ! Description:
 ! ------------
 !> Subroutine defining 3D output variables
-!------------------------------------------------------------------------------!
- SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do,  &
-                           nzt_do )
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
 
     IMPLICIT NONE
@@ -1082,4 +1085,5 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
 
+
     flag_nr  = 0
     found    = .TRUE.
@@ -1092,5 +1096,5 @@
              to_be_resorted => ti
           ELSE
-             IF ( .NOT. ALLOCATED( ti_av ) ) THEN
+             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
                 ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 ti_av = REAL( fill_value, KIND = wp )
@@ -1104,5 +1108,5 @@
              to_be_resorted => uu
           ELSE
-             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
+             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
                 ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 uu_av = REAL( fill_value, KIND = wp )
@@ -1116,5 +1120,5 @@
              to_be_resorted => vv
           ELSE
-             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
+             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
                 ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 vv_av = REAL( fill_value, KIND = wp )
@@ -1128,5 +1132,5 @@
              to_be_resorted => ww
           ELSE
-             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
+             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
                 ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 ww_av = REAL( fill_value, KIND = wp )
@@ -1140,5 +1144,5 @@
              to_be_resorted => wu
           ELSE
-             IF ( .NOT. ALLOCATED( wu_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wu_av ) )  THEN
                 ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wu_av = REAL( fill_value, KIND = wp )
@@ -1152,5 +1156,5 @@
              to_be_resorted => wv
           ELSE
-             IF ( .NOT. ALLOCATED( wv_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wv_av ) )  THEN
                 ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wv_av = REAL( fill_value, KIND = wp )
@@ -1164,5 +1168,5 @@
              to_be_resorted => wtheta
           ELSE
-             IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
                 ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wtheta_av = REAL( fill_value, KIND = wp )
@@ -1176,5 +1180,5 @@
              to_be_resorted => wq
           ELSE
-             IF ( .NOT. ALLOCATED( wq_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wq_av ) )  THEN
                 ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wq_av = REAL( fill_value, KIND = wp )
@@ -1188,5 +1192,5 @@
              to_be_resorted => wspeed
           ELSE
-             IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
                 ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wspeed_av = REAL( fill_value, KIND = wp )
@@ -1200,5 +1204,5 @@
              to_be_resorted => wdir
           ELSE
-             IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
+             IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
                 ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wdir_av = REAL( fill_value, KIND = wp )
@@ -1217,7 +1221,7 @@
           DO  j = nys, nyn
              DO  k = nzb_do, nzt_do
-                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
-                                     REAL( fill_value, KIND = wp ),            &
-                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
+                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                                    &
+                                         REAL( fill_value, KIND = wp ),                            &
+                                         BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
              ENDDO
           ENDDO
@@ -1227,9 +1231,12 @@
  END SUBROUTINE doq_output_3d
 
+
+!--------------------------------------------------------------------------------------------------!
+!
 ! Description:
 ! ------------
-!> Resorts the user-defined output quantity with indices (k,j,i) to a
-!> temporary array with indices (i,j,k) for masked data output.
-!------------------------------------------------------------------------------!
+!> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with indices
+!> (i,j,k) for masked data output.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
 
@@ -1239,4 +1246,6 @@
 
     IMPLICIT NONE
+
+    REAL(wp), PARAMETER   ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
 
     CHARACTER (LEN=*) ::  variable   !<
@@ -1257,10 +1266,8 @@
     LOGICAL      ::  resorted        !< true if array is resorted
 
-    REAL(wp),                                                                  &
-       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
-          local_pf   !<
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
-
-    REAL(wp), PARAMETER   ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
+    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf   !<
+
+    REAL(wp), DIMENSION(:,:,:), POINTER  ::  to_be_resorted  !< points to array which needs to be resorted for output
+
 
     flag_nr  = 0
@@ -1373,12 +1380,11 @@
              DO  j = 1, mask_size_l(mid,2)
                 DO  k = 1, mask_size_l(mid,3)
-                   local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k),  &
-                                                           mask_j(mid,j),  &
-                                                           mask_i(mid,i)), &
-                                            REAL( fill_value, KIND = wp ), &
-                                            BTEST( wall_flags_total_0(     &
-                                                           mask_k(mid,k),  &
-                                                           mask_j(mid,j),  &
-                                                           mask_i(mid,i)), &
+                   local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k),                          &
+                                                           mask_j(mid,j),                          &
+                                                           mask_i(mid,i)),                         &
+                                            REAL( fill_value, KIND = wp ),                         &
+                                            BTEST( wall_flags_total_0(mask_k(mid,k),               &
+                                                                      mask_j(mid,j),               &
+                                                                      mask_i(mid,i)),              &
                                                    flag_nr ) )
                 ENDDO
@@ -1395,6 +1401,5 @@
                 im = mask_i(mid,i)
                 jm = mask_j(mid,j)
-                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
-                              DIM = 1 ) - 1
+                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ), DIM=1 ) - 1
                 DO  k = 1, mask_size_l(mid,3)
                    kk = MIN( ktt+mask_k(mid,k), nzt+1 )
@@ -1415,9 +1420,9 @@
  END SUBROUTINE doq_output_mask
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Allocate required arrays
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE doq_init
 
@@ -1425,4 +1430,5 @@
 
     INTEGER(iwp) ::  ivar   !< loop index over all 2d/3d/mask output quantities
+
 
 !
@@ -1564,7 +1570,7 @@
 
 !
-!-- Save timestep number to check in time_integration if doq_calculate
-!-- has been called already, since the CALL occurs at two locations, but the calculations need to be
-!-- done only once per timestep.
+!-- Save timestep number to check in time_integration if doq_calculate has been called already,
+!-- since the CALL occurs at two locations, but the calculations need to be done only once per
+!-- timestep.
     timestep_number_at_prev_calc = current_timestep_number
 
@@ -1580,18 +1586,18 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt
-                      ti(k,j,i) = 0.25_wp * SQRT(                              &
-                        (   (   w(k,j+1,i) + w(k-1,j+1,i)                      &
-                              - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy              &
-                          - (   v(k+1,j,i) + v(k+1,j+1,i)                      &
-                              - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2     &
-                      + (   (   u(k+1,j,i) + u(k+1,j,i+1)                      &
-                              - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)          &
-                          - (   w(k,j,i+1) + w(k-1,j,i+1)                      &
-                              - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2     &
-                      + (   (   v(k,j,i+1) + v(k,j+1,i+1)                      &
-                              - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx              &
-                          - (   u(k,j+1,i) + u(k,j+1,i+1)                      &
-                              - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )  &
-                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0) )
+                      ti(k,j,i) = 0.25_wp * SQRT(                                                  &
+                                       (   (   w(k,j+1,i) + w(k-1,j+1,i)                           &
+                                             - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy                   &
+                                         - (   v(k+1,j,i) + v(k+1,j+1,i)                           &
+                                             - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2          &
+                                     + (   (   u(k+1,j,i) + u(k+1,j,i+1)                           &
+                                             - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)               &
+                                         - (   w(k,j,i+1) + w(k-1,j,i+1)                           &
+                                             - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2          &
+                                     + (   (   v(k,j,i+1) + v(k,j+1,i+1)                           &
+                                             - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx                   &
+                                         - (   u(k,j+1,i) + u(k,j+1,i+1)                           &
+                                             - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )       &
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -1603,7 +1609,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt
-                      uu(k,j,i) = u(k,j,i) * u(k,j,i)                          &
-                       * MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_total_0(k,j,i), 1) )
+                      uu(k,j,i) = u(k,j,i) * u(k,j,i)                                              &
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
                    ENDDO
                 ENDDO
@@ -1615,7 +1620,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt
-                      vv(k,j,i) = v(k,j,i) * v(k,j,i)                          &
-                       * MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_total_0(k,j,i), 2) )
+                      vv(k,j,i) = v(k,j,i) * v(k,j,i)                                              &
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
                    ENDDO
                 ENDDO
@@ -1627,7 +1631,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt-1
-                      ww(k,j,i) = w(k,j,i) * w(k,j,i)                          &
-                       * MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_total_0(k,j,i), 3) )
+                      ww(k,j,i) = w(k,j,i) * w(k,j,i)                                              &
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
                    ENDDO
                 ENDDO
@@ -1639,9 +1642,7 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt-1
-                      wu(k,j,i) = w(k,j,i)                                     &
-                       * 0.25_wp * ( u(k,j,i)   + u(k,j,i+1)                   &
-                                   + u(k+1,j,i) + u(k+1,j,i+1) )               &
-                       * MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_total_0(k,j,i), 0) )
+                      wu(k,j,i) = w(k,j,i)                                                         &
+                                  * 0.25_wp * ( u(k,j,i) + u(k,j,i+1) + u(k+1,j,i) + u(k+1,j,i+1) )&
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -1653,9 +1654,7 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt-1
-                      wv(k,j,i) = w(k,j,i)                                     &
-                       * 0.25_wp * ( v(k,j,i)   + v(k,j+1,i)                   &
-                                   + v(k+1,j,i) + v(k+1,j+1,i) )               &
-                       * MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_total_0(k,j,i), 0) )
+                      wv(k,j,i) = w(k,j,i)                                                         &
+                                  * 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k+1,j,i) + v(k+1,j+1,i) )&
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -1667,8 +1666,7 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt-1
-                      wtheta(k,j,i) = w(k,j,i)                                 &
-                       *  0.5_wp  * ( pt(k,j,i) + pt(k+1,j,i) )                &
-                       * MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_total_0(k,j,i), 3) )
+                      wtheta(k,j,i) = w(k,j,i)                                                     &
+                                     *  0.5_wp  * ( pt(k,j,i) + pt(k+1,j,i) )                      &
+                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ))
                    ENDDO
                 ENDDO
@@ -1680,7 +1678,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt-1
-                      wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) )&
-                       * MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_total_0(k,j,i), 3) )
+                      wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) )                    &
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
                    ENDDO
                 ENDDO
@@ -1690,8 +1687,7 @@
           CASE ( 'theta_2m*' )
 !
-!--          2-m potential temperature is caluclated from surface arrays. In
-!--          case the 2m level is below the first grid point, MOST is applied,
-!--          else, linear interpolation between two vertical grid levels is
-!--          applied. To access all surfaces, iterate over all horizontally-
+!--          2-m potential temperature is caluclated from surface arrays. In case the 2m level is
+!--          below the first grid point, MOST is applied, else, linear interpolation between two
+!--          vertical grid levels is applied. To access all surfaces, iterate over all horizontally-
 !--          upward facing surface types.
              surf => surf_def_h(0)
@@ -1705,9 +1701,8 @@
           CASE ( 'wspeed_10m*' )
 !
-!--          10-m wind speed is caluclated from surface arrays. In
-!--          case the 10m level is below the first grid point, MOST is applied,
-!--          else, linear interpolation between two vertical grid levels is
-!--          applied. To access all surfaces, iterate over all horizontally-
-!--          upward facing surface types.
+!--          10-m wind speed is caluclated from surface arrays. In case the 10m level is below the
+!--          first grid point, MOST is applied, else, linear interpolation between two vertical grid
+!--          levels is applied. To access all surfaces, iterate over all horizontally-upward facing
+!--          surface types.
              surf => surf_def_h(0)
              CALL calc_wind_10m
@@ -1717,12 +1712,12 @@
              CALL calc_wind_10m
 !
-!--       horizontal wind speed
+!--       Horizontal wind speed
           CASE ( 'wspeed' )
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    DO  k = nzb, nzt+1
-                      wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2             &
-                                          + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 )           &
-                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0) )
+                      wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2              &
+                                          + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 )            &
+                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ))
                    ENDDO
                 ENDDO
@@ -1730,5 +1725,5 @@
 
 !
-!--       horizontal wind direction
+!--       Horizontal wind direction
           CASE ( 'wdir' )
              DO  i = nxl, nxr
@@ -1738,6 +1733,6 @@
                       v_center(k,j,i) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
 
-                      wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) ) &
-                                  / pi * 180.0_wp + 180.0_wp
+                      wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) )                      &
+                                    / pi * 180.0_wp + 180.0_wp
                    ENDDO
                 ENDDO
@@ -1752,15 +1747,14 @@
 
 !
-!-- The following block contains subroutines to calculate diagnostic
-!-- surface quantities.
+!-- The following block contains subroutines to calculate diagnostic surface quantities.
     CONTAINS
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Calculation of 2-m potential temperature.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
        SUBROUTINE calc_pt_2m
 
-          USE surface_layer_fluxes_mod,                                        &
+          USE surface_layer_fluxes_mod,                                                            &
               ONLY:  psi_h
 
@@ -1776,17 +1770,16 @@
              k = surf%k(m)
 !
-!--          If 2-m level is below the first grid level, MOST is
-!--          used for calculation of 2-m temperature.
+!--          If 2-m level is below the first grid level, MOST is used for calculation of
+!--          2-m temperature.
              IF ( surf%z_mo(m) > 2.0_wp )  THEN
-                pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa           &
-                                * ( LOG( 2.0_wp /  surf%z0h(m) )               &
-                                  - psi_h( 2.0_wp / surf%ol(m) )               &
-                                  + psi_h( surf%z0h(m) / surf%ol(m) ) )
-!
-!--          If 2-m level is above the first grid level, 2-m temperature
-!--          is linearly interpolated between the two nearest vertical grid
-!--          levels. Note, since 2-m temperature is only computed for
-!--          horizontal upward-facing surfaces, only a vertical
-!--          interpolation is necessary.
+                pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa                               &
+                                * ( LOG( 2.0_wp / surf%z0h(m) )                                    &
+                                    - psi_h( 2.0_wp      / surf%ol(m) )                            &
+                                    + psi_h( surf%z0h(m) / surf%ol(m) ) )
+!
+!--          If 2-m level is above the first grid level, 2-m temperature is linearly interpolated
+!--          between the two nearest vertical grid levels. Note, since 2-m temperature is only
+!--          computed for horizontal upward-facing surfaces, only a vertical interpolation is
+!--          necessary.
              ELSE
 !
@@ -1798,7 +1791,7 @@
 !
 !--             kk defines the index of the first grid level >= 2m.
-                pt_2m(j,i) = pt(kk-1,j,i) +                                    &
-                              ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *            &
-                              ( pt(kk,j,i)       - pt(kk-1,j,i) ) /            &
+                pt_2m(j,i) = pt(kk-1,j,i) +                                                        &
+                              ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *                                &
+                              ( pt(kk,j,i)       - pt(kk-1,j,i) ) /                                &
                               ( zu(kk)           - zu(kk-1)     )
              ENDIF
@@ -1808,12 +1801,12 @@
        END SUBROUTINE calc_pt_2m
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Calculation of 10-m wind speed.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
        SUBROUTINE calc_wind_10m
 
-          USE surface_layer_fluxes_mod,                                        &
+          USE surface_layer_fluxes_mod,                                                            &
               ONLY:  psi_m
 
@@ -1825,4 +1818,5 @@
           REAL(wp) ::  uv_l !< wind speed at lower grid point
           REAL(wp) ::  uv_u !< wind speed at upper grid point
+
 
           DO  m = 1, surf%ns
@@ -1832,17 +1826,16 @@
              k = surf%k(m)
 !
-!--          If 10-m level is below the first grid level, MOST is
-!--          used for calculation of 10-m temperature.
+!--          If 10-m level is below the first grid level, MOST is used for calculation of 10-m
+!--          temperature.
              IF ( surf%z_mo(m) > 10.0_wp )  THEN
-                uv_10m(j,i) = surf%us(m) / kappa                               &
-                          * ( LOG( 10.0_wp /  surf%z0(m) )                     &
-                              - psi_m( 10.0_wp    / surf%ol(m) )               &
-                              + psi_m( surf%z0(m) / surf%ol(m) ) )
-!
-!--          If 10-m level is above the first grid level, 10-m wind speed
-!--          is linearly interpolated between the two nearest vertical grid
-!--          levels. Note, since 10-m temperature is only computed for
-!--          horizontal upward-facing surfaces, only a vertical
-!--          interpolation is necessary.
+                uv_10m(j,i) = surf%us(m) / kappa                                                   &
+                              * ( LOG( 10.0_wp /  surf%z0(m) )                                     &
+                                  - psi_m( 10.0_wp    / surf%ol(m) )                               &
+                                  + psi_m( surf%z0(m) / surf%ol(m) ) )
+!
+!--          If 10-m level is above the first grid level, 10-m wind speed is linearly interpolated
+!--          between the two nearest vertical grid levels. Note, since 10-m temperature is only
+!--          computed for horizontal upward-facing surfaces, only a vertical interpolation is
+!--          necessary.
              ELSE
 !
@@ -1854,12 +1847,12 @@
 !
 !--             kk defines the index of the first grid level >= 10m.
-                uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2   &
+                uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2                       &
                            + ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
 
-                uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i)   + u(kk,j,i+1)   ) )**2   &
+                uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i)   + u(kk,j,i+1)   ) )**2                       &
                            + ( 0.5_wp * ( v(kk,j,i)   + v(kk,j+1,i)   ) )**2 )
 
-                uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) *        &
-                                     ( uv_u              - uv_l     ) /        &
+                uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) *                            &
+                                     ( uv_u              - uv_l     ) /                            &
                                      ( zu(kk)            - zu(kk-1) )
 
@@ -1873,19 +1866,18 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Preparation of the diagnostic output, counting of the module-specific
-!> output quantities and gathering of the output names.
-!------------------------------------------------------------------------------!
+!> Preparation of the diagnostic output, counting of the module-specific output quantities and
+!> gathering of the output names.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE doq_prepare
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  do2d, do3d, domask, masks
 
     IMPLICIT NONE
 
-    CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d_var = ' '  !<
-                                                          !< label array for 2d output quantities
+    CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d_var = ' '  !< label array for 2d output quantities
 
     INTEGER(iwp) ::  av         !< index defining type of output, av=0 instantaneous, av=1 averaged
@@ -1895,4 +1887,5 @@
     INTEGER(iwp) ::  mid          !< masked output running index
 
+
     prepared_diagnostic_output_quantities = .FALSE.
 
@@ -1907,6 +1900,5 @@
 !
 !--    Gather 2d output quantity names.
-!--    Check for double occurrence of output quantity, e.g. by _xy,
-!--    _yz, _xz.
+!--    Check for double occurrence of output quantity, e.g. by _xy, _yz, _xz.
        DO  WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
           IF ( .NOT.  ANY( do_all == do2d_var(av,ivar) ) )  THEN
@@ -1930,6 +1922,5 @@
        ivar = 1
 !
-!--    Gather masked output quantity names. Also check for double output
-!--    e.g. by different masks.
+!--    Gather masked output quantity names. Also check for double output e.g. by different masks.
        DO  mid = 1, masks
           DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
@@ -1950,12 +1941,11 @@
  END SUBROUTINE doq_prepare
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine reads local (subdomain) restart data
-!------------------------------------------------------------------------------!
- SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,                             &
-                               nxr_on_file, nynf, nync, nyn_on_file, nysf,                         &
-                               nysc, nys_on_file, tmp_2d, tmp_3d, found )
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,    &
+                               nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, found )
 
     USE control_parameters
@@ -2002,5 +1992,5 @@
           pt_2m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =                                     &
                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-    
+
        CASE ( 'ti_av' )
           IF ( .NOT. ALLOCATED( ti_av ) )  THEN
@@ -2108,9 +2098,9 @@
  END SUBROUTINE doq_rrd_local_ftn
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Read module-specific local restart data arrays (MPI-IO).
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE doq_rrd_local_mpi
 
@@ -2198,13 +2188,13 @@
  END SUBROUTINE doq_rrd_local_mpi
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine writes local (subdomain) restart data
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE doq_wrd_local
 
-
     IMPLICIT NONE
+
 
     IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
Index: palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_s.f90	(revision 4582)
+++ palm/trunk/SOURCE/diffusion_s.f90	(revision 4583)
@@ -1,20 +1,19 @@
 !> @file diffusion_s.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:
@@ -25,19 +24,22 @@
 ! -----------------
 ! $Id$
-! Introduction of wall_flags_total_0, which currently sets bits based on static
-! topography information used in wall_flags_static_0
-! 
+! file re-formatted to follow the PALM coding standard
+!
+! 4360 2020-01-07 11:25:50Z suehring
+! Introduction of wall_flags_total_0, which currently sets bits based on static topography
+! information used in wall_flags_static_0
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_total_0
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 3927 2019-04-23 13:24:29Z raasch
 ! pointer attribute removed from scalar 3d-array for performance reasons
-! 
+!
 ! 3761 2019-02-25 15:31:42Z raasch
 ! unused variables removed
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! nopointer option removed
@@ -50,7 +52,7 @@
 ! ------------
 !> Diffusion term of scalar quantities (temperature and water content)
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  MODULE diffusion_s_mod
- 
+
 
     PRIVATE
@@ -65,37 +67,35 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for all grid points
-!------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_s( s, s_flux_def_h_up,    s_flux_def_h_down,          &
-                               s_flux_t,                                       &
-                               s_flux_lsm_h_up,    s_flux_usm_h_up,            &
-                               s_flux_def_v_north, s_flux_def_v_south,         &
-                               s_flux_def_v_east,  s_flux_def_v_west,          &
-                               s_flux_lsm_v_north, s_flux_lsm_v_south,         &
-                               s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
-                               s_flux_usm_v_north, s_flux_usm_v_south,         &
+!--------------------------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_s( s, s_flux_def_h_up,    s_flux_def_h_down,                              &
+                               s_flux_t,                                                           &
+                               s_flux_lsm_h_up,    s_flux_usm_h_up,                                &
+                               s_flux_def_v_north, s_flux_def_v_south,                             &
+                               s_flux_def_v_east,  s_flux_def_v_west,                              &
+                               s_flux_lsm_v_north, s_flux_lsm_v_south,                             &
+                               s_flux_lsm_v_east,  s_flux_lsm_v_west,                              &
+                               s_flux_usm_v_north, s_flux_usm_v_south,                             &
                                s_flux_usm_v_east,  s_flux_usm_v_west )
 
-       USE arrays_3d,                                                          &
+       USE arrays_3d,                                                                              &
            ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
-       
-       USE control_parameters,                                                 & 
+
+       USE control_parameters,                                                                     &
            ONLY: use_surface_fluxes, use_top_fluxes
-       
-       USE grid_variables,                                                     &
+
+       USE grid_variables,                                                                         &
            ONLY:  ddx, ddx2, ddy, ddy2
-       
-       USE indices,                                                            &
-           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,        &
-                  wall_flags_total_0
-       
+
+       USE indices,                                                                                &
+           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0
+
        USE kinds
 
-       USE surface_mod,                                                        &
-           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
-                   surf_usm_v 
+       USE surface_mod,                                                                            &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
        IMPLICIT NONE
@@ -109,26 +109,26 @@
 
        REAL(wp) ::  flag              !< flag to mask topography grid points
-       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
-       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point 
+       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface
+       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
        REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
-       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point 
+       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
+       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface
        REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
-       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface  
-
+
+       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
        REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
        REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
-       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
        REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
-       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
-       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces 
        REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical natural-type surfaces
        REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
        REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
-       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical natural-type surfaces
        REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical natural-type surfaces
        REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
-       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
@@ -164,8 +164,8 @@
 !
 !--             Predetermine flag to mask topography and wall-bounded grid points
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 
-!
-!--             Predetermine flag to mask wall-bounded grid points, equivalent to
-!--             former s_outer array
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--             Predetermine flag to mask wall-bounded grid points, equivalent to former s_outer
+!--             array
                 mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
                 mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
@@ -173,26 +173,25 @@
                 mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
 
-                tend(k,j,i) = tend(k,j,i)                                      &
-                                          + 0.5_wp * (                         &
-                        mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )               &
-                                   * ( s(k,j,i+1) - s(k,j,i)   )               &
-                      - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )               &
-                                   * ( s(k,j,i)   - s(k,j,i-1) )               &
-                                                     ) * ddx2 * flag           &
-                                          + 0.5_wp * (                         &
-                        mask_north * ( kh(k,j,i) + kh(k,j+1,i) )               &
-                                   * ( s(k,j+1,i) - s(k,j,i)   )               &
-                      - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )               &
-                                   * ( s(k,j,i)   - s(k,j-1,i) )               &
-                                                     ) * ddy2 * flag
-             ENDDO
-
-!
-!--          Apply prescribed horizontal wall heatflux where necessary. First,
-!--          determine start and end index for respective (j,i)-index. Please
-!--          note, in the flat case following loop will not be entered, as
-!--          surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
-!--          so far. 
-!--          First, for default-type surfaces 
+                tend(k,j,i) = tend(k,j,i)                                                          &
+                                                  + 0.5_wp * (                                     &
+                                mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )                           &
+                                           * ( s(k,j,i+1) - s(k,j,i)   )                           &
+                              - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )                           &
+                                           * ( s(k,j,i)   - s(k,j,i-1) )                           &
+                                                             ) * ddx2 * flag                       &
+                                                  + 0.5_wp * (                                     &
+                                mask_north * ( kh(k,j,i) + kh(k,j+1,i) )                           &
+                                           * ( s(k,j+1,i) - s(k,j,i)   )                           &
+                              - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )                           &
+                                           * ( s(k,j,i)   - s(k,j-1,i) )                           &
+                                                             ) * ddy2 * flag
+             ENDDO
+
+!
+!--          Apply prescribed horizontal wall heatflux where necessary. First, determine start and
+!--          end index for respective (j,i)-index. Please note, in the flat case following loop will
+!--          not be entered, as surf_s=1 and surf_e=0. Furtermore, note, no vertical natural
+!--          surfaces so far.
+!--          First, for default-type surfaces.
 !--          North-facing vertical default-type surfaces
              surf_s = surf_def_v(0)%start_index(j,i)
@@ -294,34 +293,28 @@
 
 !
-!--          Compute vertical diffusion. In case that surface fluxes have been
-!--          prescribed or computed at bottom and/or top, index k starts/ends at
-!--          nzb+2 or nzt-1, respectively. Model top is also mask if top flux
-!--          is given.
+!--          Compute vertical diffusion. In case that surface fluxes have been prescribed or
+!--          computed at bottom and/or top, index k starts/ends at nzb+2 or nzt-1, respectively.
+!--          Model top is also mask if top flux is given.
              DO  k = nzb+1, nzt
 !
-!--             Determine flags to mask topography below and above. Flag 0 is 
-!--             used to mask topography in general, and flag 8 implies 
-!--             information about use_surface_fluxes. Flag 9 is used to control 
-!--             flux at model top. 
-                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 
-                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
-                              MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k+1,j,i), 9 ) ) 
-                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k,j,i), 0 ) )
-
-                tend(k,j,i) = tend(k,j,i)                                      &
-                                       + 0.5_wp * (                            &
-                                      ( kh(k,j,i) + kh(k+1,j,i) ) *            &
-                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
-                                                            * rho_air_zw(k)    &
-                                                            * mask_top         &
-                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
-                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                                            * rho_air_zw(k-1)  &
-                                                            * mask_bottom      &
-                                                  ) * ddzw(k) * drho_air(k)    &
+!--             Determine flags to mask topography below and above. Flag 0 is used to mask
+!--             topography in general, and flag 8 implies information about use_surface_fluxes.
+!--             Flag 9 is used to control flux at model top.
+                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
+                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
+                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
+                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+
+                tend(k,j,i) = tend(k,j,i)                                                          &
+                                       + 0.5_wp * (                                                &
+                                      ( kh(k,j,i) + kh(k+1,j,i) ) *                                &
+                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)                      &
+                                                            * rho_air_zw(k)                        &
+                                                            * mask_top                             &
+                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *                                &
+                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)                        &
+                                                            * rho_air_zw(k-1)                      &
+                                                            * mask_bottom                          &
+                                                  ) * ddzw(k) * drho_air(k)                        &
                                                               * flag
              ENDDO
@@ -331,5 +324,5 @@
              IF ( use_surface_fluxes )  THEN
 !
-!--             Default-type surfaces, upward-facing               
+!--             Default-type surfaces, upward-facing
                 surf_s = surf_def_h(0)%start_index(j,i)
                 surf_e = surf_def_h(0)%end_index(j,i)
@@ -337,10 +330,9 @@
 
                    k   = surf_def_h(0)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)              &
-                                       * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) * ddzw(k) * drho_air(k)
 
                 ENDDO
 !
-!--             Default-type surfaces, downward-facing               
+!--             Default-type surfaces, downward-facing
                 surf_s = surf_def_h(1)%start_index(j,i)
                 surf_e = surf_def_h(1)%end_index(j,i)
@@ -348,10 +340,9 @@
 
                    k   = surf_def_h(1)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)            &
-                                       * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) * ddzw(k) * drho_air(k)
 
                 ENDDO
 !
-!--             Natural-type surfaces, upward-facing  
+!--             Natural-type surfaces, upward-facing
                 surf_s = surf_lsm_h%start_index(j,i)
                 surf_e = surf_lsm_h%end_index(j,i)
@@ -359,10 +350,9 @@
 
                    k   = surf_lsm_h%k(m)
-                   tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m)              &
-                                       * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) * ddzw(k) * drho_air(k)
 
                 ENDDO
 !
-!--             Urban-type surfaces, upward-facing     
+!--             Urban-type surfaces, upward-facing
                 surf_s = surf_usm_h%start_index(j,i)
                 surf_e = surf_usm_h%end_index(j,i)
@@ -370,6 +360,5 @@
 
                    k   = surf_usm_h%k(m)
-                   tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m)              &
-                                       * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) * ddzw(k) * drho_air(k)
 
                 ENDDO
@@ -384,6 +373,5 @@
 
                    k   = surf_def_h(2)%k(m)
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                           + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
                 ENDDO
              ENDIF
@@ -394,37 +382,36 @@
     END SUBROUTINE diffusion_s
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_s_ij( i, j, s,                                        &
-                               s_flux_def_h_up,    s_flux_def_h_down,          &
-                               s_flux_t,                                       &
-                               s_flux_lsm_h_up,    s_flux_usm_h_up,            &
-                               s_flux_def_v_north, s_flux_def_v_south,         &
-                               s_flux_def_v_east,  s_flux_def_v_west,          &
-                               s_flux_lsm_v_north, s_flux_lsm_v_south,         &
-                               s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
-                               s_flux_usm_v_north, s_flux_usm_v_south,         &
-                               s_flux_usm_v_east,  s_flux_usm_v_west )       
-
-       USE arrays_3d,                                                          &
+!--------------------------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_s_ij( i, j, s,                                                            &
+                               s_flux_def_h_up,    s_flux_def_h_down,                              &
+                               s_flux_t,                                                           &
+                               s_flux_lsm_h_up,    s_flux_usm_h_up,                                &
+                               s_flux_def_v_north, s_flux_def_v_south,                             &
+                               s_flux_def_v_east,  s_flux_def_v_west,                              &
+                               s_flux_lsm_v_north, s_flux_lsm_v_south,                             &
+                               s_flux_lsm_v_east,  s_flux_lsm_v_west,                              &
+                               s_flux_usm_v_north, s_flux_usm_v_south,                             &
+                               s_flux_usm_v_east,  s_flux_usm_v_west )
+
+       USE arrays_3d,                                                                              &
            ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
-           
-       USE control_parameters,                                                 & 
+
+       USE control_parameters,                                                                     &
            ONLY: use_surface_fluxes, use_top_fluxes
-       
-       USE grid_variables,                                                     &
+
+       USE grid_variables,                                                                         &
            ONLY:  ddx, ddx2, ddy, ddy2
-       
-       USE indices,                                                            &
+
+       USE indices,                                                                                &
            ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_total_0
-       
+
        USE kinds
 
-       USE surface_mod,                                                        &
-           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
-                   surf_usm_v 
+       USE surface_mod,                                                                            &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
        IMPLICIT NONE
@@ -438,26 +425,26 @@
 
        REAL(wp) ::  flag              !< flag to mask topography grid points
-       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
-       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point 
+       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface
+       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
        REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
-       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point 
+       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
+       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface
        REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
-       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface  
-
+
+       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
        REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
        REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
-       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
        REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
-       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
-       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces 
        REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
-       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
-       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
        REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
@@ -465,4 +452,5 @@
        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !< treated scalar
 
+
 !
 !--    Compute horizontal diffusion
@@ -470,8 +458,7 @@
 !
 !--       Predetermine flag to mask topography and wall-bounded grid points
-          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 
-!
-!--       Predetermine flag to mask wall-bounded grid points, equivalent to
-!--       former s_outer array
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--       Predetermine flag to mask wall-bounded grid points, equivalent to former s_outer array
           mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
           mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
@@ -479,29 +466,27 @@
           mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
 !
-!--       Finally, determine flag to mask both topography itself as well
-!--       as wall-bounded grid points, which will be treated further below
-
-          tend(k,j,i) = tend(k,j,i)                                            &
-                                          + 0.5_wp * (                         &
-                            mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )           &
-                                       * ( s(k,j,i+1) - s(k,j,i)   )           &
-                          - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )           &
-                                       * ( s(k,j,i)   - s(k,j,i-1) )           &
-                                                     ) * ddx2 * flag           &
-                                          + 0.5_wp * (                         &
-                            mask_north * ( kh(k,j,i) + kh(k,j+1,i) )           &
-                                       * ( s(k,j+1,i) - s(k,j,i)   )           &
-                          - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )           &
-                                       * ( s(k,j,i)  - s(k,j-1,i)  )           &
+!--       Finally, determine flag to mask both topography itself as well as wall-bounded grid
+!--       points, which will be treated further below
+
+          tend(k,j,i) = tend(k,j,i)                                                                &
+                                          + 0.5_wp * (                                             &
+                            mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )                               &
+                                       * ( s(k,j,i+1) - s(k,j,i)   )                               &
+                          - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )                               &
+                                       * ( s(k,j,i)   - s(k,j,i-1) )                               &
+                                                     ) * ddx2 * flag                               &
+                                          + 0.5_wp * (                                             &
+                            mask_north * ( kh(k,j,i) + kh(k,j+1,i) )                               &
+                                       * ( s(k,j+1,i) - s(k,j,i)   )                               &
+                          - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )                               &
+                                       * ( s(k,j,i)  - s(k,j-1,i)  )                               &
                                                      ) * ddy2 * flag
        ENDDO
 
 !
-!--    Apply prescribed horizontal wall heatflux where necessary. First,
-!--    determine start and end index for respective (j,i)-index. Please
-!--    note, in the flat case following loops will not be entered, as
-!--    surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
-!--    so far. 
-!--    First, for default-type surfaces
+!--    Apply prescribed horizontal wall heatflux where necessary. First, determine start and end
+!--    index for respective (j,i)-index. Please note, in the flat case following loops will not be
+!--    entered, as surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces so far.
+!--    First, for default-type surfaces.
 !--    North-facing vertical default-type surfaces
        surf_s = surf_def_v(0)%start_index(j,i)
@@ -604,34 +589,28 @@
 
 !
-!--    Compute vertical diffusion. In case that surface fluxes have been
-!--    prescribed or computed at bottom and/or top, index k starts/ends at
-!--    nzb+2 or nzt-1, respectively. Model top is also mask if top flux 
-!--    is given. 
+!--    Compute vertical diffusion. In case that surface fluxes have been prescribed or computed at
+!--    bottom and/or top, index k starts/ends at nzb+2 or nzt-1, respectively. Model top is also
+!--    mask if top flux is given.
        DO  k = nzb+1, nzt
 !
-!--       Determine flags to mask topography below and above. Flag 0 is 
-!--       used to mask topography in general, and flag 8 implies 
-!--       information about use_surface_fluxes. Flag 9 is used to control 
-!--       flux at model top.   
-          mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 
-          mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k+1,j,i), 8 ) )  *    &
-                        MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
-          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
-
-          tend(k,j,i) = tend(k,j,i)                                            &
-                                       + 0.5_wp * (                            &
-                                      ( kh(k,j,i) + kh(k+1,j,i) ) *            &
-                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
-                                                            * rho_air_zw(k)    &
-                                                            * mask_top         &
-                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
-                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                                            * rho_air_zw(k-1)  &
-                                                            * mask_bottom      &
-                                                  ) * ddzw(k) * drho_air(k)    &
+!--       Determine flags to mask topography below and above. Flag 0 is used to mask topography in
+!--       general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
+!--       control flux at model top.
+          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
+          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) )  *        &
+                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
+          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+
+          tend(k,j,i) = tend(k,j,i)                                                                &
+                                       + 0.5_wp * (                                                &
+                                      ( kh(k,j,i) + kh(k+1,j,i) ) *                                &
+                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)                      &
+                                                            * rho_air_zw(k)                        &
+                                                            * mask_top                             &
+                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *                                &
+                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)                        &
+                                                            * rho_air_zw(k-1)                      &
+                                                            * mask_bottom                          &
+                                                  ) * ddzw(k) * drho_air(k)                        &
                                                               * flag
        ENDDO
@@ -649,6 +628,5 @@
              k   = surf_def_h(0)%k(m)
 
-             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)                    &
-                                       * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) * ddzw(k) * drho_air(k)
           ENDDO
 !
@@ -660,6 +638,5 @@
              k   = surf_def_h(1)%k(m)
 
-             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)                  &
-                                       * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) * ddzw(k) * drho_air(k)
           ENDDO
 !
@@ -670,6 +647,5 @@
              k   = surf_lsm_h%k(m)
 
-             tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m)                    &
-                                       * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) * ddzw(k) * drho_air(k)
           ENDDO
 !
@@ -680,6 +656,5 @@
              k   = surf_usm_h%k(m)
 
-             tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m)                    &
-                                       * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) * ddzw(k) * drho_air(k)
           ENDDO
        ENDIF
@@ -692,6 +667,5 @@
 
              k   = surf_def_h(2)%k(m)
-             tend(k,j,i) = tend(k,j,i)                                         &
-                           + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
           ENDDO
        ENDIF
Index: palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_u.f90	(revision 4582)
+++ palm/trunk/SOURCE/diffusion_u.f90	(revision 4583)
@@ -1,20 +1,19 @@
 !> @file diffusion_u.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:
@@ -25,13 +24,16 @@
 ! -----------------
 ! $Id$
-! Introduction of wall_flags_total_0, which currently sets bits based on static
-! topography information used in wall_flags_static_0
-! 
+! file re-formatted to follow the PALM coding standard
+!
+! 4360 2020-01-07 11:25:50Z suehring
+! Introduction of wall_flags_total_0, which currently sets bits based on static topography
+! information used in wall_flags_static_0
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! OpenACC port for SPEC
@@ -44,9 +46,9 @@
 ! ------------
 !> Diffusion term of the u-component
-!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization
-!>       and slows down the speed on NEC about 5-10%
-!------------------------------------------------------------------------------!
+!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization and slows down the
+!        speed on NEC about 5-10%
+!--------------------------------------------------------------------------------------------------!
  MODULE diffusion_u_mod
- 
+
 
     PRIVATE
@@ -61,29 +63,27 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for all grid points
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
     SUBROUTINE diffusion_u
 
-       USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
-       
-       USE control_parameters,                                                 &
-           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
-                  use_top_fluxes
-       
-       USE grid_variables,                                                     &
+       USE arrays_3d,                                                                              &
+           ONLY:  ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
+
+       USE control_parameters,                                                                     &
+           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
+
+       USE grid_variables,                                                                         &
            ONLY:  ddx, ddx2, ddy
-       
-       USE indices,                                                            &
+
+       USE indices,                                                                                &
            ONLY:  nxlu, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
-     
+
        USE kinds
 
-       USE surface_mod,                                                        &
-           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
-                   surf_usm_v
+       USE surface_mod,                                                                            &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
        IMPLICIT NONE
@@ -102,8 +102,8 @@
        REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
        REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
-       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface       
-       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point 
-       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point 
-       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface 
+       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
+       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
+       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
+       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
 
 
@@ -125,37 +125,34 @@
              DO  k = nzb+1, nzt
 !
-!--             Predetermine flag to mask topography and wall-bounded grid points. 
-!--             It is sufficient to masked only north- and south-facing surfaces, which
-!--             need special treatment for the u-component. 
-                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) ) 
+!--             Predetermine flag to mask topography and wall-bounded grid points.
+!--             It is sufficient to masked only north- and south-facing surfaces, which need special
+!--             treatment for the u-component.
+                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) )
                 mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
                 mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmyp = 0.25_wp *                                               &
-                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-                kmym = 0.25_wp *                                               &
-                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
-
-                tend(k,j,i) = tend(k,j,i)                                      &
-                        + 2.0_wp * (                                           &
-                                  km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
-                                - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
-                                   ) * ddx2 * flag                             &
-                        +          ( mask_north * (                            &
-                            kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
-                          + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
-                                                  )                            &
-                                   - mask_south * (                            &
-                            kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
-                          + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
-                                                  )                            &
-                                   ) * ddy  * flag                             
+                kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+                kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
+
+                tend(k,j,i) = tend(k,j,i)                                                          &
+                              + 2.0_wp * (                                                         &
+                                        km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )                  &
+                                      - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )                  &
+                                         ) * ddx2 * flag                                           &
+                              +          ( mask_north * (                                          &
+                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy                       &
+                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx                       &
+                                                        )                                          &
+                                         - mask_south * (                                          &
+                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy                           &
+                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx                           &
+                                                        )                                          &
+                                         ) * ddy  * flag
              ENDDO
 !
-!--          Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1)
-!--          surfaces. Note, in the the flat case, loops won't be entered as 
-!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
-!--          so far.           
+!--          Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces.
+!--          Note, in the the flat case, loops won't be entered as start_index > end_index.
+!--          Furtermore, note, no vertical natural surfaces so far.
 !--          Default-type surfaces
              DO  l = 0, 1
@@ -164,7 +161,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_def_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &                    
-                                    surf_def_v(l)%mom_flux_uv(m) * ddy
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
+                ENDDO
              ENDDO
 !
@@ -175,7 +171,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_lsm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &                    
-                                    surf_lsm_v(l)%mom_flux_uv(m) * ddy
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
+                ENDDO
              ENDDO
 !
@@ -186,55 +181,45 @@
                 DO  m = surf_s, surf_e
                    k           = surf_usm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &                    
-                                    surf_usm_v(l)%mom_flux_uv(m) * ddy
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
+                ENDDO
              ENDDO
 
 !
-!--          Compute vertical diffusion. In case of simulating a surface layer,
-!--          respective grid diffusive fluxes are masked (flag 8) within this 
-!--          loop, and added further below, else, simple gradient approach is
-!--          applied. Model top is also mask if top-momentum flux is given. 
+!--          Compute vertical diffusion. In case of simulating a surface layer, respective grid
+!--          diffusive fluxes are masked (flag 8) within this loop, and added further below, else,
+!--          simple gradient approach is applied. Model top is also mask if top-momentum flux is
+!--          given.
              DO  k = nzb+1, nzt
 !
-!--             Determine flags to mask topography below and above. Flag 1 is 
-!--             used to mask topography in general, and flag 8 implies 
-!--             information about use_surface_fluxes. Flag 9 is used to control 
-!--             momentum flux at model top.  
-                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 
-                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
-                              MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k+1,j,i), 9 ) ) 
-                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k,j,i), 1 ) ) 
+!--             Determine flags to mask topography below and above. Flag 1 is used to mask
+!--             topography in general, and flag 8 implies information about use_surface_fluxes.
+!--             Flag 9 is used to control momentum flux at model top.
+                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
+                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
+                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
+                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25_wp *                                               &
-                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-                kmzm = 0.25_wp *                                               &
-                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
-
-                tend(k,j,i) = tend(k,j,i)                                      &
-                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
-                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
-                                   ) * rho_air_zw(k)   * mask_top              &
-                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
-                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
-                                   ) * rho_air_zw(k-1) * mask_bottom           &
-                          ) * ddzw(k) * drho_air(k) * flag
+                kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
+                kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+
+                tend(k,j,i) = tend(k,j,i)                                                          &
+                              + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)                 &
+                                         + ( w(k,j,i)   - w(k,j,i-1) ) * ddx                       &
+                                         ) * rho_air_zw(k)   * mask_top                            &
+                                - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)                 &
+                                         + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx                     &
+                                         ) * rho_air_zw(k-1) * mask_bottom                         &
+                                ) * ddzw(k) * drho_air(k) * flag
              ENDDO
 
 !
-!--          Vertical diffusion at the first grid point above the surface,
-!--          if the momentum flux at the bottom is given by the Prandtl law or
-!--          if it is prescribed by the user.
-!--          Difference quotient of the momentum flux is not formed over half
-!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
-!--          with other (LES) models showed that the values of the momentum
-!--          flux becomes too large in this case.
-!--          The term containing w(k-1,..) (see above equation) is removed here
-!--          because the vertical velocity is assumed to be zero at the surface.
+!--          Vertical diffusion at the first grid point above the surface, if the momentum flux at
+!--          the bottom is given by the Prandtl law or if it is prescribed by the user.
+!--          Difference quotient of the momentum flux is not formed over half of the grid spacing
+!--          (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
+!--          values of the momentum flux becomes too large in this case.
+!--          The term containing w(k-1,..) (see above equation) is removed here because the vertical
+!--          velocity is assumed to be zero at the surface.
              IF ( use_surface_fluxes )  THEN
 !
@@ -246,7 +231,6 @@
                    k   = surf_def_h(0)%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                        + ( - ( - surf_def_h(0)%usws(m) )                      &
-                          ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
                 ENDDO
 !
@@ -258,7 +242,6 @@
                    k   = surf_def_h(1)%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                        + ( - surf_def_h(1)%usws(m)                            &
-                          ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
                 ENDDO
 !
@@ -270,7 +253,6 @@
                    k   = surf_lsm_h%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                        + ( - ( - surf_lsm_h%usws(m) )                         &
-                          ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - ( - surf_lsm_h%usws(m) ) ) * ddzw(k) * drho_air(k)
                 ENDDO
 !
@@ -282,7 +264,6 @@
                    k   = surf_usm_h%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                       + ( - ( - surf_usm_h%usws(m) )                          &
-                         ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - ( - surf_usm_h%usws(m) ) ) * ddzw(k) * drho_air(k)
                 ENDDO
 
@@ -297,6 +278,6 @@
                    k   = surf_def_h(2)%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                        + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
                 ENDDO
              ENDIF
@@ -308,29 +289,27 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for grid point i,j
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
     SUBROUTINE diffusion_u_ij( i, j )
 
-       USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
-       
-       USE control_parameters,                                                 &
-           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
-                  use_top_fluxes
-       
-       USE grid_variables,                                                     &
+       USE arrays_3d,                                                                              &
+           ONLY:  ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
+
+       USE control_parameters,                                                                     &
+           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
+
+       USE grid_variables,                                                                         &
            ONLY:  ddx, ddx2, ddy
-       
-       USE indices,                                                            &
+
+       USE indices,                                                                                &
            ONLY:  nzb, nzt, wall_flags_total_0
-      
+
        USE kinds
 
-       USE surface_mod,                                                        &
-           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
-                   surf_usm_v
+       USE surface_mod,                                                                            &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
        IMPLICIT NONE
@@ -349,16 +328,16 @@
        REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
        REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
-       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface       
-       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point 
-       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point 
-       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface 
-! 
+       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
+       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
+       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
+       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
+!
 !--    Compute horizontal diffusion
        DO  k = nzb+1, nzt
 !
-!--       Predetermine flag to mask topography and wall-bounded grid points. 
-!--       It is sufficient to masked only north- and south-facing surfaces, which
-!--       need special treatment for the u-component. 
-          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) ) 
+!--       Predetermine flag to mask topography and wall-bounded grid points.
+!--       It is sufficient to masked only north- and south-facing surfaces, which need special
+!--       treatment for the u-component.
+          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) )
           mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
           mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
@@ -368,24 +347,23 @@
           kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
-          tend(k,j,i) = tend(k,j,i)                                            &
-                       + 2.0_wp * (                                            &
-                                 km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )     &
-                               - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )     &
-                                   ) * ddx2 * flag                             &
-                                 + (                                           &
-                  mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i)     ) * ddy    &
-                                      + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx    &
-                                      )                                        &
-                - mask_south * kmym * ( ( u(k,j,i)   - u(k,j-1,i)   ) * ddy    &
-                                      + ( v(k,j,i)   - v(k,j,i-1)   ) * ddx    &
-                                      )                                        &
-                                   ) * ddy  * flag
-       ENDDO
-
-!
-!--    Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1)
-!--    surfaces. Note, in the the flat case, loops won't be entered as 
-!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
-!--    so far.           
+          tend(k,j,i) = tend(k,j,i)                                                                &
+                                + 2.0_wp * (                                                       &
+                                         km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )                 &
+                                       - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )                 &
+                                           ) * ddx2 * flag                                         &
+                                         + (                                                       &
+                          mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i)     ) * ddy                &
+                                              + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx                &
+                                              )                                                    &
+                        - mask_south * kmym * ( ( u(k,j,i)   - u(k,j-1,i)   ) * ddy                &
+                                              + ( v(k,j,i)   - v(k,j,i-1)   ) * ddx                &
+                                              )                                                    &
+                                           ) * ddy  * flag
+       ENDDO
+
+!
+!--    Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces. Note, in
+!--    the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
+!--    vertical natural surfaces so far.
 !--    Default-type surfaces
        DO  l = 0, 1
@@ -395,5 +373,5 @@
              k           = surf_def_v(l)%k(m)
              tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
-          ENDDO   
+          ENDDO
        ENDDO
 !
@@ -405,5 +383,5 @@
              k           = surf_lsm_v(l)%k(m)
              tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
-          ENDDO   
+          ENDDO
        ENDDO
 !
@@ -415,25 +393,19 @@
              k           = surf_usm_v(l)%k(m)
              tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
-          ENDDO   
-       ENDDO
-!
-!--    Compute vertical diffusion. In case of simulating a surface layer,
-!--    respective grid diffusive fluxes are masked (flag 8) within this 
-!--    loop, and added further below, else, simple gradient approach is
-!--    applied. Model top is also mask if top-momentum flux is given.
+          ENDDO
+       ENDDO
+!
+!--    Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
+!--    fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
+!--    approach is applied. Model top is also mask if top-momentum flux is given.
        DO  k = nzb+1, nzt
 !
-!--       Determine flags to mask topography below and above. Flag 1 is 
-!--       used to mask topography in general, and flag 8 implies 
-!--       information about use_surface_fluxes. Flag 9 is used to control 
-!--       momentum flux at model top. 
-          mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 
-          mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *     &
-                        MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
-          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k,j,i), 1 ) ) 
+!--       Determine flags to mask topography below and above. Flag 1 is used to mask topography in
+!--       general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
+!--       control momentum flux at model top.
+          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
+          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *         &
+                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
+          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
@@ -441,22 +413,20 @@
           kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
-          tend(k,j,i) = tend(k,j,i)                                            &
-                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
-                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
-                                   ) * rho_air_zw(k)   * mask_top              &
-                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
-                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
-                                   ) * rho_air_zw(k-1) * mask_bottom           &
+          tend(k,j,i) = tend(k,j,i)                                                                &
+                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)                       &
+                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx                             &
+                                   ) * rho_air_zw(k)   * mask_top                                  &
+                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)                       &
+                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx                           &
+                                   ) * rho_air_zw(k-1) * mask_bottom                               &
                           ) * ddzw(k) * drho_air(k) * flag
        ENDDO
 
 !
-!--    Vertical diffusion at the first surface grid points, if the
-!--    momentum flux at the bottom is given by the Prandtl law or if it is
-!--    prescribed by the user.
-!--    Difference quotient of the momentum flux is not formed over half of
-!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with 
-!--    other (LES) models showed that the values of the momentum flux becomes
-!--    too large in this case.
+!--    Vertical diffusion at the first surface grid points, if the momentum flux at the bottom is
+!--    given by the Prandtl law or if it is prescribed by the user.
+!--    Difference quotient of the momentum flux is not formed over half of the grid spacing
+!--    (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
+!--    of the momentum flux becomes too large in this case.
        IF ( use_surface_fluxes )  THEN
 !
@@ -468,7 +438,5 @@
              k   = surf_def_h(0)%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                        + ( - ( - surf_def_h(0)%usws(m) )                      &
-                          ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
           ENDDO
 !
@@ -480,7 +448,5 @@
              k   = surf_def_h(1)%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                        + ( - surf_def_h(1)%usws(m)                            &
-                          ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
           ENDDO
 !
@@ -492,7 +458,5 @@
              k   = surf_lsm_h%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                        + ( - ( - surf_lsm_h%usws(m) )                         &
-                          ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h%usws(m) ) ) * ddzw(k) * drho_air(k)
           ENDDO
 !
@@ -504,7 +468,5 @@
              k   = surf_usm_h%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                       + ( - ( - surf_usm_h%usws(m) )                          &
-                         ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h%usws(m) ) ) * ddzw(k) * drho_air(k)
           ENDDO
 
@@ -519,6 +481,5 @@
              k   = surf_def_h(2)%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                           + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
           ENDDO
        ENDIF
Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 4582)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 4583)
@@ -1,20 +1,19 @@
 !> @file diffusion_v.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:
@@ -25,13 +24,16 @@
 ! -----------------
 ! $Id$
-! Introduction of wall_flags_total_0, which currently sets bits based on static
-! topography information used in wall_flags_static_0
-! 
+! file re-formatted to follow the PALM coding standard
+!
+! 4360 2020-01-07 11:25:50Z suehring
+! Introduction of wall_flags_total_0, which currently sets bits based on static topography
+! information used in wall_flags_static_0
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! OpenACC port for SPEC
@@ -44,7 +46,7 @@
 ! ------------
 !> Diffusion term of the v-component
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  MODULE diffusion_v_mod
- 
+
 
     PRIVATE
@@ -59,29 +61,27 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for all grid points
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
     SUBROUTINE diffusion_v
 
-       USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
-       
-       USE control_parameters,                                                 &
-           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
-                  use_top_fluxes
-       
-       USE grid_variables,                                                     &
+       USE arrays_3d,                                                                              &
+           ONLY:  ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w 
+
+       USE control_parameters,                                                                     &
+           ONLY:  constant_top_momentumflux, use_surface_fluxes,  use_top_fluxes
+
+       USE grid_variables,                                                                         &
            ONLY:  ddx, ddy, ddy2
-       
-       USE indices,                                                            &
+
+       USE indices,                                                                                &
            ONLY:  nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_total_0
-       
+
        USE kinds
 
-       USE surface_mod,                                                        &
-           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
-                   surf_usm_v
+       USE surface_mod,                                                                            &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
        IMPLICIT NONE
@@ -100,8 +100,8 @@
        REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto yv-zw grid
        REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto yv-zw grid
-       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface  
-       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point 
-       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point 
-       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface      
+       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
+       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
+       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
+       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
 
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
@@ -122,8 +122,8 @@
 
 !
-!--             Predetermine flag to mask topography and wall-bounded grid points. 
-!--             It is sufficient to masked only east- and west-facing surfaces, which
-!--             need special treatment for the v-component. 
-                flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) ) 
+!--             Predetermine flag to mask topography and wall-bounded grid points.
+!--             It is sufficient to masked only east- and west-facing surfaces, which need special
+!--             treatment for the v-component.
+                flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) )
                 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
                 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
@@ -133,26 +133,25 @@
                 kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
-                tend(k,j,i) = tend(k,j,i) +    (                             &
-                          mask_east * kmxp * (                               &
-                                 ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
-                               + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
-                                             )                               &
-                        - mask_west * kmxm * (                               &
-                                 ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
-                               + ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
-                                             )                               &
-                                               ) * ddx  * flag               &
-                                    + 2.0_wp * (                             &
-                                  km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )  &
-                                - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )  &
-                                               ) * ddy2 * flag
+                tend(k,j,i) = tend(k,j,i) +          (                                             &
+                                mask_east * kmxp * (                                               &
+                                       ( v(k,j,i+1) - v(k,j,i)     ) * ddx                         &
+                                     + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy                         &
+                                                   )                                               &
+                              - mask_west * kmxm * (                                               &
+                                       ( v(k,j,i) - v(k,j,i-1) ) * ddx                             &
+                                     + ( u(k,j,i) - u(k,j-1,i) ) * ddy                             &
+                                                   )                                               &
+                                                     ) * ddx  * flag                               &
+                                          + 2.0_wp * (                                             &
+                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )                  &
+                                      - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )                  &
+                                                     ) * ddy2 * flag
 
              ENDDO
 
 !
-!--          Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
-!--          surfaces. Note, in the the flat case, loops won't be entered as 
-!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
-!--          so far.           
+!--          Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note,
+!--          in the the flat case, loops won't be entered as start_index > end_index. Furtermore,
+!--          note, no vertical natural surfaces so far.
 !--          Default-type surfaces
              DO  l = 2, 3
@@ -161,7 +160,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_def_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                    surf_def_v(l)%mom_flux_uv(m) * ddx
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
+                ENDDO
              ENDDO
 !
@@ -172,7 +170,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_lsm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                    surf_lsm_v(l)%mom_flux_uv(m) * ddx
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
+                ENDDO
              ENDDO
 !
@@ -183,52 +180,42 @@
                 DO  m = surf_s, surf_e
                    k           = surf_usm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                    surf_usm_v(l)%mom_flux_uv(m) * ddx
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
+                ENDDO
              ENDDO
 !
-!--          Compute vertical diffusion. In case of simulating a surface layer,
-!--          respective grid diffusive fluxes are masked (flag 10) within this 
-!--          loop, and added further below, else, simple gradient approach is
-!--          applied. Model top is also mask if top-momentum flux is given. 
+!--          Compute vertical diffusion. In case of simulating a surface layer, respective grid
+!--          diffusive fluxes are masked (flag 10) within this loop, and added further below, else,
+!--          simple gradient approach is applied. Model top is also mask if top-momentum flux is
+!--          given.
              DO  k = nzb+1, nzt
 !
-!--             Determine flags to mask topography below and above. Flag 2 is 
-!--             used to mask topography in general, while flag 8 implies also
-!--             information about use_surface_fluxes. Flag 9 is used to control 
-!--             momentum flux at model top.  
-                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 
-                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
-                              MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k+1,j,i), 9 ) ) 
-                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
-                                 BTEST( wall_flags_total_0(k,j,i), 2 ) ) 
+!--             Determine flags to mask topography below and above. Flag 2 is used to mask
+!--             topography in general, while flag 8 implies also information about
+!--             use_surface_fluxes. Flag 9 is used to control momentum flux at model top.
+                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
+                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
+                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
+                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25_wp * &
-                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                kmzm = 0.25_wp * &
-                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
-
-                tend(k,j,i) = tend(k,j,i)                                      &
-                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
-                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
-                      &            ) * rho_air_zw(k)   * mask_top              &
-                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
-                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
-                      &            ) * rho_air_zw(k-1) * mask_bottom           &
-                      &   ) * ddzw(k) * drho_air(k) * flag
+                kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                tend(k,j,i) = tend(k,j,i)                                                          &
+                              & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)                 &
+                              &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy                       &
+                              &            ) * rho_air_zw(k)   * mask_top                          &
+                              &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)               &
+                              &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy                   &
+                              &            ) * rho_air_zw(k-1) * mask_bottom                       &
+                              &   ) * ddzw(k) * drho_air(k) * flag
              ENDDO
 
 !
-!--          Vertical diffusion at the first grid point above the surface,
-!--          if the momentum flux at the bottom is given by the Prandtl law
-!--          or if it is prescribed by the user.
-!--          Difference quotient of the momentum flux is not formed over
-!--          half of the grid spacing (2.0*ddzw(k)) any more, since the
-!--          comparison with other (LES) models showed that the values of
-!--          the momentum flux becomes too large in this case.
+!--          Vertical diffusion at the first grid point above the surface, if the momentum flux at
+!--          the bottom is given by the Prandtl law or if it is prescribed by the user.
+!--          Difference quotient of the momentum flux is not formed over half of the grid spacing
+!--          (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
+!--          values of the momentum flux becomes too large in this case.
              IF ( use_surface_fluxes )  THEN
 !
@@ -239,7 +226,6 @@
                    k   = surf_def_h(0)%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                        + ( - ( - surf_def_h(0)%vsws(m) )                      &
-                          ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
                 ENDDO
 !
@@ -250,7 +236,6 @@
                    k   = surf_def_h(1)%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                        + ( - surf_def_h(1)%vsws(m)                            &
-                          ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
                 ENDDO
 !
@@ -261,7 +246,6 @@
                    k   = surf_lsm_h%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                        + ( - ( - surf_lsm_h%vsws(m) )                         &
-                          ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - ( - surf_lsm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
 
                 ENDDO
@@ -273,7 +257,6 @@
                    k   = surf_usm_h%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                        + ( - ( - surf_usm_h%vsws(m) )                         &
-                          ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - ( - surf_usm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
 
                 ENDDO
@@ -288,6 +271,6 @@
                    k   = surf_def_h(2)%k(m)
 
-                   tend(k,j,i) = tend(k,j,i)                                   &
-                           + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
+                   tend(k,j,i) = tend(k,j,i)                                                       &
+                                 + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
                 ENDDO
              ENDIF
@@ -299,29 +282,27 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for grid point i,j
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
     SUBROUTINE diffusion_v_ij( i, j )
 
-       USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
-       
-       USE control_parameters,                                                 &
-           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
-                  use_top_fluxes
-       
-       USE grid_variables,                                                     &
+       USE arrays_3d,                                                                              &
+           ONLY:  ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
+
+       USE control_parameters,                                                                     &
+           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
+
+       USE grid_variables,                                                                         &
            ONLY:  ddx, ddy, ddy2
-       
-       USE indices,                                                            &
+
+       USE indices,                                                                                &
            ONLY:  nzb, nzt, wall_flags_total_0
-       
+
        USE kinds
 
-       USE surface_mod,                                                        &
-           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
-                   surf_usm_v
+       USE surface_mod,                                                                            &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
        IMPLICIT NONE
@@ -341,7 +322,7 @@
        REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
        REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
-       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface  
-       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point 
-       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point 
+       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
+       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
+       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
        REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
 
@@ -350,8 +331,8 @@
        DO  k = nzb+1, nzt
 !
-!--       Predetermine flag to mask topography and wall-bounded grid points. 
-!--       It is sufficient to masked only east- and west-facing surfaces, which
-!--       need special treatment for the v-component. 
-          flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) ) 
+!--       Predetermine flag to mask topography and wall-bounded grid points.
+!--       It is sufficient to masked only east- and west-facing surfaces, which need special
+!--       treatment for the v-component.
+          flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) )
           mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
           mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
@@ -361,25 +342,24 @@
           kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
-          tend(k,j,i) = tend(k,j,i) +          (                             &
-                          mask_east * kmxp * (                               &
-                                 ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
-                               + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
-                                             )                               &
-                        - mask_west * kmxm * (                               &
-                                 ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
-                               + ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
-                                             )                               &
-                                               ) * ddx  * flag               &
-                                    + 2.0_wp * (                             &
-                                  km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )  &
-                                - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )  &
+          tend(k,j,i) = tend(k,j,i) +          (                                                   &
+                          mask_east * kmxp * (                                                     &
+                                 ( v(k,j,i+1) - v(k,j,i)     ) * ddx                               &
+                               + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy                               &
+                                             )                                                     &
+                        - mask_west * kmxm * (                                                     &
+                                 ( v(k,j,i) - v(k,j,i-1) ) * ddx                                   &
+                               + ( u(k,j,i) - u(k,j-1,i) ) * ddy                                   &
+                                             )                                                     &
+                                               ) * ddx  * flag                                     &
+                                    + 2.0_wp * (                                                   &
+                                  km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )                        &
+                                - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )                        &
                                                ) * ddy2 * flag
        ENDDO
 
 !
-!--    Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
-!--    surfaces. Note, in the the flat case, loops won't be entered as 
-!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
-!--    so far.           
+!--    Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note, in the
+!--    the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
+!--    vertical natural surfaces so far.
 !--    Default-type surfaces
        DO  l = 2, 3
@@ -389,5 +369,5 @@
              k           = surf_def_v(l)%k(m)
              tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
-          ENDDO   
+          ENDDO
        ENDDO
 !
@@ -399,5 +379,5 @@
              k           = surf_lsm_v(l)%k(m)
              tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
-          ENDDO   
+          ENDDO
        ENDDO
 !
@@ -409,25 +389,19 @@
              k           = surf_usm_v(l)%k(m)
              tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
-          ENDDO   
-       ENDDO
-!
-!--    Compute vertical diffusion. In case of simulating a surface layer,
-!--    respective grid diffusive fluxes are masked (flag 8) within this 
-!--    loop, and added further below, else, simple gradient approach is
-!--    applied. Model top is also mask if top-momentum flux is given.
+          ENDDO
+       ENDDO
+!
+!--    Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
+!--    fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
+!--    approach is applied. Model top is also mask if top-momentum flux is given.
        DO  k = nzb+1, nzt
 !
-!--       Determine flags to mask topography below and above. Flag 2 is 
-!--       used to mask topography in general, while flag 10 implies also
-!--       information about use_surface_fluxes. Flag 9 is used to control 
-!--       momentum flux at model top.  
-          mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 
-          mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *     &
-                        MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k+1,j,i), 9 ) ) 
-          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_total_0(k,j,i), 2 ) )
+!--       Determine flags to mask topography below and above. Flag 2 is used to mask topography in
+!--       general, while flag 10 implies also information about use_surface_fluxes. Flag 9 is used
+!--       to control momentum flux at model top.
+          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
+          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *         &
+                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
+          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
@@ -435,22 +409,20 @@
           kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
-          tend(k,j,i) = tend(k,j,i)                                            &
-                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
-                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
-                      &            ) * rho_air_zw(k)   * mask_top              &
-                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
-                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
-                      &            ) * rho_air_zw(k-1) * mask_bottom           &
-                      &   ) * ddzw(k) * drho_air(k) * flag
-       ENDDO
-
-!
-!--    Vertical diffusion at the first grid point above the surface, if the
-!--    momentum flux at the bottom is given by the Prandtl law or if it is
-!--    prescribed by the user.
-!--    Difference quotient of the momentum flux is not formed over half of
-!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with 
-!--    other (LES) models showed that the values of the momentum flux becomes
-!--    too large in this case.
+          tend(k,j,i) = tend(k,j,i)                                                                &
+                        & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)                       &
+                        &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy                             &
+                        &            ) * rho_air_zw(k)   * mask_top                                &
+                        &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)                     &
+                        &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy                         &
+                        &            ) * rho_air_zw(k-1) * mask_bottom                             &
+                        &   ) * ddzw(k) * drho_air(k) * flag
+       ENDDO
+
+!
+!--    Vertical diffusion at the first grid point above the surface, if the momentum flux at the
+!--    bottom is given by the Prandtl law or if it is prescribed by the user.
+!--    Difference quotient of the momentum flux is not formed over half of the grid spacing
+!--    (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
+!--    of the momentum flux becomes too large in this case.
        IF ( use_surface_fluxes )  THEN
 !
@@ -461,7 +433,5 @@
              k   = surf_def_h(0)%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                        + ( - ( - surf_def_h(0)%vsws(m) )                      &
-                          ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
           ENDDO
 !
@@ -472,7 +442,5 @@
              k   = surf_def_h(1)%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                        + ( - surf_def_h(1)%vsws(m)                            &
-                          ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
           ENDDO
 !
@@ -483,7 +451,5 @@
              k   = surf_lsm_h%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                        + ( - ( - surf_lsm_h%vsws(m) )                         &
-                          ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
 
           ENDDO
@@ -495,7 +461,5 @@
              k   = surf_usm_h%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                        + ( - ( - surf_usm_h%vsws(m) )                         &
-                          ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
 
           ENDDO
@@ -510,6 +474,5 @@
              k   = surf_def_h(2)%k(m)
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                           + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
+             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
           ENDDO
        ENDIF
Index: palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_w.f90	(revision 4582)
+++ palm/trunk/SOURCE/diffusion_w.f90	(revision 4583)
@@ -1,20 +1,19 @@
 !> @file diffusion_w.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:
@@ -25,13 +24,16 @@
 ! -----------------
 ! $Id$
-! Introduction of wall_flags_total_0, which currently sets bits based on static
-! topography information used in wall_flags_static_0
-! 
+! file re-formatted to follow the PALM coding standard
+!
+! 4360 2020-01-07 11:25:50Z suehring
+! Introduction of wall_flags_total_0, which currently sets bits based on static topography
+! information used in wall_flags_static_0
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! OpenACC port for SPEC
@@ -44,7 +46,7 @@
 ! ------------
 !> Diffusion term of the w-component
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  MODULE diffusion_w_mod
- 
+
 
     PRIVATE
@@ -59,23 +61,23 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for all grid points
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
     SUBROUTINE diffusion_w
 
-       USE arrays_3d,                                                          &          
-           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
-           
-       USE grid_variables,                                                     &     
+       USE arrays_3d,                                                                              &
+           ONLY :  ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
+
+       USE grid_variables,                                                                         &
            ONLY :  ddx, ddy
-           
-       USE indices,                                                            &            
+
+       USE indices,                                                                                &
            ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
-           
+
        USE kinds
 
-       USE surface_mod,                                                        &
+       USE surface_mod,                                                                            &
            ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
 
@@ -89,5 +91,5 @@
        INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
        INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
-       
+
        REAL(wp) ::  flag              !< flag to mask topography grid points
        REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
@@ -95,8 +97,8 @@
        REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
        REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
+       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
+       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
+       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
        REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
-       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point 
-       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point 
-       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point 
 
 
@@ -116,58 +118,52 @@
              DO  k = nzb+1, nzt-1
 !
-!--             Predetermine flag to mask topography and wall-bounded grid points. 
-                flag       = MERGE( 1.0_wp, 0.0_wp,                            &
-                                    BTEST( wall_flags_total_0(k,j,i),   3 ) ) 
-                mask_east  = MERGE( 1.0_wp, 0.0_wp,                            &
-                                    BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
-                mask_west  = MERGE( 1.0_wp, 0.0_wp,                            &
-                                    BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
-                mask_south = MERGE( 1.0_wp, 0.0_wp,                            &
-                                    BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
-                mask_north = MERGE( 1.0_wp, 0.0_wp,                            &
-                                    BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
+!--             Predetermine flag to mask topography and wall-bounded grid points.
+                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) )
+                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
+                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
+                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
+                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +               &
+                kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +                                   &
                                    km(k+1,j,i) + km(k+1,j,i+1) )
-                kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +               &
+                kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +                                   &
                                    km(k+1,j,i) + km(k+1,j,i-1) )
-                kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
+                kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +                                   &
                                    km(k,j+1,i) + km(k+1,j+1,i) )
-                kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
+                kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +                                   &
                                    km(k,j-1,i) + km(k+1,j-1,i) )
 
-                tend(k,j,i) = tend(k,j,i)                                      &
-                       + ( mask_east *  kmxp * (                               &
-                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
-                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
-                                               )                               &
-                         - mask_west * kmxm *  (                               &
-                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
-                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
-                                               )                               &
-                         ) * ddx                                 * flag        &
-                       + ( mask_north * kmyp * (                               &
-                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
-                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
-                                               )                               &
-                         - mask_south * kmym * (                               &
-                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
-                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
-                                               )                               &
-                         ) * ddy                                 * flag        &
-                       + 2.0_wp * (                                            &
-                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1) &
-                                     * rho_air(k+1)                            &
-                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
-                                     * rho_air(k)                              &
-                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
-             ENDDO
-
-!
-!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
-!--          surfaces. Note, in the the flat case, loops won't be entered as 
-!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
-!--          so far.           
+                tend(k,j,i) = tend(k,j,i)                                                          &
+                              + ( mask_east *  kmxp * (                                            &
+                                          ( w(k,j,i+1)   - w(k,j,i)   ) * ddx                      &
+                                        + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)                &
+                                                      )                                            &
+                                - mask_west * kmxm *  (                                            &
+                                          ( w(k,j,i)     - w(k,j,i-1) ) * ddx                      &
+                                        + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)                &
+                                                      )                                            &
+                                ) * ddx                                 * flag                     &
+                              + ( mask_north * kmyp * (                                            &
+                                          ( w(k,j+1,i)   - w(k,j,i)   ) * ddy                      &
+                                        + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)                &
+                                                      )                                            &
+                                - mask_south * kmym * (                                            &
+                                          ( w(k,j,i)     - w(k,j-1,i) ) * ddy                      &
+                                        + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)                &
+                                                      )                                            &
+                                ) * ddy                                 * flag                     &
+                              + 2.0_wp * (                                                         &
+                                km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1)              &
+                                            * rho_air(k+1)                                         &
+                              - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)                &
+                                            * rho_air(k)                                           &
+                                         ) * ddzu(k+1) * drho_air_zw(k) * flag
+             ENDDO
+
+!
+!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces.
+!--          Note, in the the flat case, loops won't be entered as start_index > end_index.
+!--          Furtermore, note, no vertical natural surfaces so far.
 !--          Default-type surfaces
              DO  l = 0, 1
@@ -176,7 +172,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_def_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                     surf_def_v(l)%mom_flux_w(m) * ddy
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
+                ENDDO
              ENDDO
 !
@@ -187,7 +182,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_lsm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                     surf_lsm_v(l)%mom_flux_w(m) * ddy
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
+                ENDDO
              ENDDO
 !
@@ -198,11 +192,9 @@
                 DO  m = surf_s, surf_e
                    k           = surf_usm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                     surf_usm_v(l)%mom_flux_w(m) * ddy
-                ENDDO   
-             ENDDO
-!
-!--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
-!--          surface.
+                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
+                ENDDO
+             ENDDO
+!
+!--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surface.
 !--          Default-type surfaces
              DO  l = 2, 3
@@ -211,7 +203,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_def_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                     surf_def_v(l)%mom_flux_w(m) * ddx
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
+                ENDDO
              ENDDO
 !
@@ -222,7 +213,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_lsm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                     surf_lsm_v(l)%mom_flux_w(m) * ddx
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
+                ENDDO
              ENDDO
 !
@@ -233,7 +223,6 @@
                 DO  m = surf_s, surf_e
                    k           = surf_usm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) +                                 &
-                                     surf_usm_v(l)%mom_flux_w(m) * ddx
-                ENDDO   
+                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
+                ENDDO
              ENDDO
 
@@ -244,23 +233,23 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for grid point i,j
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
     SUBROUTINE diffusion_w_ij( i, j )
 
-       USE arrays_3d,                                                          &          
-           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
-           
-       USE grid_variables,                                                     &     
+       USE arrays_3d,                                                                              &
+           ONLY :  ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
+
+       USE grid_variables,                                                                         &
            ONLY :  ddx, ddy
-           
-       USE indices,                                                            &            
+
+       USE indices,                                                                                &
            ONLY :  nzb, nzt, wall_flags_total_0
-           
+
        USE kinds
 
-       USE surface_mod,                                                        &
+       USE surface_mod,                                                                            &
            ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
 
@@ -275,5 +264,5 @@
        INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
        INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
-       
+
        REAL(wp) ::  flag              !< flag to mask topography grid points
        REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
@@ -281,14 +270,14 @@
        REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
        REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
+       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
+       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
+       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
        REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
-       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point 
-       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point 
-       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point 
 
 
        DO  k = nzb+1, nzt-1
 !
-!--       Predetermine flag to mask topography and wall-bounded grid points. 
-          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) ) 
+!--       Predetermine flag to mask topography and wall-bounded grid points.
+          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) )
           mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
           mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
@@ -302,35 +291,34 @@
           kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
-          tend(k,j,i) = tend(k,j,i)                                            &
-                       + ( mask_east *  kmxp * (                               &
-                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
-                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
-                                               )                               &
-                         - mask_west * kmxm *  (                               &
-                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
-                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
-                                               )                               &
-                         ) * ddx                                 * flag        &
-                       + ( mask_north * kmyp * (                               &
-                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
-                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
-                                               )                               &
-                         - mask_south * kmym * (                               &
-                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
-                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
-                                               )                               &
-                         ) * ddy                                 * flag        &
-                       + 2.0_wp * (                                            &
-                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)   &
-                                     * rho_air(k+1)                            &
-                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
-                                     * rho_air(k)                              &
-                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
-       ENDDO
-!
-!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
-!--    surfaces. Note, in the the flat case, loops won't be entered as 
-!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
-!--    so far.           
+          tend(k,j,i) = tend(k,j,i)                                                                &
+                        + ( mask_east *  kmxp * (                                                  &
+                                    ( w(k,j,i+1)   - w(k,j,i)   ) * ddx                            &
+                                  + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)                      &
+                                                )                                                  &
+                          - mask_west * kmxm *  (                                                  &
+                                    ( w(k,j,i)     - w(k,j,i-1) ) * ddx                            &
+                                  + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)                      &
+                                                )                                                  &
+                          ) * ddx                                 * flag                           &
+                        + ( mask_north * kmyp * (                                                  &
+                                    ( w(k,j+1,i)   - w(k,j,i)   ) * ddy                            &
+                                  + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)                      &
+                                                )                                                  &
+                          - mask_south * kmym * (                                                  &
+                                    ( w(k,j,i)     - w(k,j-1,i) ) * ddy                            &
+                                  + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)                      &
+                                                )                                                  &
+                          ) * ddy                                 * flag                           &
+                        + 2.0_wp * (                                                               &
+                          km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)                      &
+                                      * rho_air(k+1)                                               &
+                        - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)                      &
+                                      * rho_air(k)                                                 &
+                                   ) * ddzu(k+1) * drho_air_zw(k) * flag
+       ENDDO
+!
+!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces. Note, in
+!--    the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
+!--    vertical natural surfaces so far.
 !--    Default-type surfaces
        DO  l = 0, 1
@@ -339,7 +327,6 @@
           DO  m = surf_s, surf_e
              k           = surf_def_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) +                                       &
-                                     surf_def_v(l)%mom_flux_w(m) * ddy
-          ENDDO   
+             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
+          ENDDO
        ENDDO
 !
@@ -350,7 +337,6 @@
           DO  m = surf_s, surf_e
              k           = surf_lsm_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) +                                       &
-                                     surf_lsm_v(l)%mom_flux_w(m) * ddy
-          ENDDO   
+             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
+          ENDDO
        ENDDO
 !
@@ -361,11 +347,9 @@
           DO  m = surf_s, surf_e
              k           = surf_usm_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) +                                       &
-                                     surf_usm_v(l)%mom_flux_w(m) * ddy
-          ENDDO   
-       ENDDO
-!
-!--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
-!--    surfaces. 
+             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
+          ENDDO
+       ENDDO
+!
+!--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surfaces.
 !--    Default-type surfaces
        DO  l = 2, 3
@@ -374,7 +358,6 @@
           DO  m = surf_s, surf_e
              k           = surf_def_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) +                                       &
-                                     surf_def_v(l)%mom_flux_w(m) * ddx
-          ENDDO   
+             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
+          ENDDO
        ENDDO
 !
@@ -385,7 +368,6 @@
           DO  m = surf_s, surf_e
              k           = surf_lsm_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) +                                       &
-                                     surf_lsm_v(l)%mom_flux_w(m) * ddx
-          ENDDO   
+             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
+          ENDDO
        ENDDO
 !
@@ -396,7 +378,6 @@
           DO  m = surf_s, surf_e
              k           = surf_usm_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) +                                       &
-                                     surf_usm_v(l)%mom_flux_w(m) * ddx
-          ENDDO   
+             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
+          ENDDO
        ENDDO
 
Index: palm/trunk/SOURCE/disturb_field.f90
===================================================================
--- palm/trunk/SOURCE/disturb_field.f90	(revision 4582)
+++ palm/trunk/SOURCE/disturb_field.f90	(revision 4583)
@@ -1,20 +1,19 @@
 !> @file disturb_field.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:
@@ -25,19 +24,22 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4457 2020-03-11 14:20:43Z raasch
 ! use statement for exchange horiz added
-! 
+!
 ! 4360 2020-01-07 11:25:50Z suehring
-! Introduction of wall_flags_total_0, which currently sets bits based on static
-! topography information used in wall_flags_static_0
-! 
+! Introduction of wall_flags_total_0, which currently sets bits based on static topography
+! information used in wall_flags_static_0
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4237 2019-09-25 11:33:42Z knoop
 ! Added missing OpenMP directives
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 3849 2019-04-01 16:35:16Z knoop
 ! Corrected "Former revisions" section
@@ -50,36 +52,33 @@
 ! ------------
 !> Imposing a random perturbation on a 3D-array.
-!> On parallel computers, the random number generator is as well called for all
-!> gridpoints of the total domain to ensure, regardless of the number of PEs
-!> used, that the elements of the array have the same values in the same
-!> order in every case. The perturbation range is steered by dist_range.
-!------------------------------------------------------------------------------!
+!> On parallel computers, the random number generator is as well called for all gridpoints of the
+!> total domain to ensure, regardless of the number of PEs used, that the elements of the array have
+!> the same values in the same order in every case. The perturbation range is steered by dist_range.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE disturb_field( var_char, dist1, field )
- 
-
-    USE control_parameters,                                                    &
-        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range,             &
-               disturbance_amplitude, disturbance_created,                     &
-               disturbance_level_ind_b, disturbance_level_ind_t, iran,         &
+
+
+    USE control_parameters,                                                                        &
+        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range, disturbance_amplitude,          &
+               disturbance_created, disturbance_level_ind_b, disturbance_level_ind_t, iran,        &
                random_generator, topography
-                
-    USE cpulog,                                                                &
+
+    USE cpulog,                                                                                    &
         ONLY:  cpu_log, log_point
-        
-    USE exchange_horiz_mod,                                                    &
+
+    USE exchange_horiz_mod,                                                                        &
         ONLY:  exchange_horiz
 
-    USE indices,                                                               &
-        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
-               nzt, wall_flags_total_0
-        
+    USE indices,                                                                                   &
+        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, nzt,                &
+               wall_flags_total_0
+
     USE kinds
-    
-    USE random_function_mod,                                                   &
+
+    USE random_function_mod,                                                                       &
         ONLY: random_function
-        
-    USE random_generator_parallel,                                             &
-        ONLY:  random_number_parallel, random_seed_parallel, random_dummy,     &
-               seq_random_array
+
+    USE random_generator_parallel,                                                                 &
+        ONLY:  random_number_parallel, random_seed_parallel, random_dummy, seq_random_array
 
     IMPLICIT NONE
@@ -87,5 +86,5 @@
     CHARACTER (LEN = *) ::  var_char !< flag to distinguish betwenn u- and v-component
 
-    INTEGER(iwp) ::  flag_nr !< number of respective flag for u- or v-grid 
+    INTEGER(iwp) ::  flag_nr !< number of respective flag for u- or v-grid
     INTEGER(iwp) ::  i       !< index variable
     INTEGER(iwp) ::  j       !< index variable
@@ -93,8 +92,8 @@
 
     REAL(wp) ::  randomnumber  !<
-    
+
     REAL(wp) ::  dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
     REAL(wp) ::  field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
-    
+
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist2  !<
 
@@ -105,6 +104,5 @@
     flag_nr = MERGE( 20, 21, TRIM(var_char) == 'u' )
 !
-!-- Create an additional temporary array and initialize the arrays needed
-!-- to store the disturbance
+!-- Create an additional temporary array and initialize the arrays needed to store the disturbance
     ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 !$ACC DATA CREATE(dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
@@ -123,9 +121,6 @@
           DO  j = dist_nys(dist_range), dist_nyn(dist_range)
              DO  k = disturbance_level_ind_b, disturbance_level_ind_t
-                randomnumber = 3.0_wp * disturbance_amplitude *                &
-                               ( random_function( iran ) - 0.5_wp )
-                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
-                     nyn >= j )                                                &
-                THEN
+                randomnumber = 3.0_wp * disturbance_amplitude * ( random_function( iran ) - 0.5_wp )
+                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND. nyn >= j )  THEN
                    dist1(k,j,i) = randomnumber
                 ENDIF
@@ -139,9 +134,6 @@
              DO  k = disturbance_level_ind_b, disturbance_level_ind_t
                 CALL random_number_parallel( random_dummy )
-                randomnumber = 3.0_wp * disturbance_amplitude *                &
-                               ( random_dummy - 0.5_wp )
-                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
-                     nyn >= j )                                                &
-                THEN
+                randomnumber = 3.0_wp * disturbance_amplitude * ( random_dummy - 0.5_wp )
+                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  nyn >= j )  THEN
                    dist1(k,j,i) = randomnumber
                 ENDIF
@@ -155,8 +147,6 @@
              DO  k = disturbance_level_ind_b, disturbance_level_ind_t
                 CALL RANDOM_NUMBER( randomnumber )
-                randomnumber = 3.0_wp * disturbance_amplitude *                &
-                                ( randomnumber - 0.5_wp )
-                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j )   &
-                THEN
+                randomnumber = 3.0_wp * disturbance_amplitude * ( randomnumber - 0.5_wp )
+                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  nyn >= j )  THEN
                    dist1(k,j,i) = randomnumber
                 ENDIF
@@ -177,8 +167,6 @@
 !
 !-- Applying the Shuman filter in order to smooth the perturbations.
-!-- Neighboured grid points in all three directions are used for the
-!-- filter operation.
-!-- Loop has been splitted to make runs reproducible on HLRN systems using
-!-- compiler option -O3
+!-- Neighboured grid points in all three directions are used for the filter operation.
+!-- Loop has been splitted to make runs reproducible on HLRN systems using compiler option -O3
      !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) PRESENT(dist1, dist2)
      !$OMP PARALLEL DO PRIVATE(i, j, k)
@@ -186,12 +174,12 @@
         DO  j = nys, nyn
           DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
-             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                  &
-                            + dist1(k,j+1,i) + dist1(k+1,j,i)                  &
+             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                                      &
+                            + dist1(k,j+1,i) + dist1(k+1,j,i)                                      &
                             ) / 12.0_wp
           ENDDO
           DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
-              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
-                            + 6.0_wp * dist1(k,j,i)                            &
-                            ) / 12.0_wp
+              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)                      &
+                            + 6.0_wp * dist1(k,j,i)                                                &
+                                            ) / 12.0_wp
           ENDDO
         ENDDO
@@ -208,7 +196,7 @@
        DO  j = nys, nyn
           DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
-             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
-                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
-                            + 6.0_wp * dist2(k,j,i)                            &
+             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i)                     &
+                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i)                     &
+                            + 6.0_wp * dist2(k,j,i)                                                &
                             ) / 12.0_wp
           ENDDO
@@ -219,13 +207,12 @@
 
 !
-!-- Remove perturbations below topography (including one gridpoint above it
-!-- in order to allow for larger timesteps at the beginning of the simulation
-!-- (diffusion criterion))
+!-- Remove perturbations below topography (including one gridpoint above it in order to allow for
+!-- larger timesteps at the beginning of the simulation (diffusion criterion))
     IF ( TRIM( topography ) /= 'flat' )  THEN
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
              DO  k = nzb, nzb_max
-                dist1(k,j,i) = MERGE( dist1(k,j,i), 0.0_wp,                    &
-                                  BTEST( wall_flags_total_0(k,j,i), flag_nr )  &
+                dist1(k,j,i) = MERGE( dist1(k,j,i), 0.0_wp,                                        &
+                                  BTEST( wall_flags_total_0(k,j,i), flag_nr )                      &
                                     )
              ENDDO
