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