Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4560)
+++ 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
