Index: palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_w.f90	(revision 4901)
+++ palm/trunk/SOURCE/diffusion_w.f90	(revision 1)
@@ -1,41 +1,52 @@
-!> @file diffusion_w.f90
-!--------------------------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
-! Public License as published by the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
-! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
-! Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with PALM. If not, see
-! <http://www.gnu.org/licenses/>.
-!
-! Copyright 1997-2021 Leibniz Universitaet Hannover
-!--------------------------------------------------------------------------------------------------!
-!
-! Current revisions:
+ MODULE diffusion_w_mod
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
 ! -----------------
 ! 
-! 
+!
 ! Former revisions:
 ! -----------------
-! $Id$
-! 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
+! $Log: diffusion_w.f90,v $
+! Revision 1.12  2006/02/23 10:38:03  raasch
+! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
+! +z0 in argument list
+! WARNING: loops containing the MAX function are still not properly vectorized!
+!
+! Revision 1.11  2005/03/26 20:11:05  raasch
+! Additional damping layer at the outflow in case of non-cyclic lateral
+! boundaries, additional arguments km_damp_x, km_damp_y
+!
+! Revision 1.10  2004/01/30 10:22:32  raasch
+! Scalar lower k index nzb replaced by 2d-array nzb_2d
+!
+! Revision 1.9  2003/03/12 16:26:49  raasch
+! Full code replaced in the call for all gridpoints instead of calling the
+! _ij version (required by NEC, because otherwise no vectorization)
+!
+! Revision 1.8  2002/06/11 12:54:16  raasch
+! Former subroutine changed to a module which allows to be called for all grid
+! points of a single vertical column with index i,j or for all grid points by
+! using function overloading.
+!
+! Revision 1.7  2001/03/30 07:12:35  raasch
+! Translation of remaining German identifiers (variables, subroutines, etc.)
+!
+! Revision 1.6  2001/01/22 06:20:12  raasch
+! Module test_variables removed
+!
+! Revision 1.5  2000/07/03 12:58:40  raasch
+! dummy arguments, whose corresponding actual arguments are pointers,
+! are now also defined as pointers,
+! all comments translated into English
+!
+! Revision 1.4  1998/07/06 12:12:24  raasch
+! + USE test_variables
+!
+! Revision 1.3  1997/09/12 07:24:58  raasch
+! Leerzeilen mussten entfernt werden
+!
+! Revision 1.2  1997/09/12 06:44:00  raasch
+! HP-Compiler erfordert & am Beginn von Fortsetzungszeilen
 !
 ! Revision 1.1  1997/09/12 06:24:11  raasch
@@ -45,8 +56,6 @@
 ! Description:
 ! ------------
-!> Diffusion term of the w-component
-!--------------------------------------------------------------------------------------------------!
- MODULE diffusion_w_mod
-
+! Diffusion term of the w-component
+!------------------------------------------------------------------------------!
 
     PRIVATE
@@ -61,169 +70,145 @@
 
 
-!--------------------------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Call for all grid points
-!--------------------------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_w
-
-       USE arrays_3d,                                                                              &
-           ONLY :  ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
-
-       USE grid_variables,                                                                         &
-           ONLY :  ddx, ddy
-
-       USE indices,                                                                                &
-           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
-
-       USE kinds
-
-       USE surface_mod,                                                                            &
-           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
+!------------------------------------------------------------------------------!
+! Call for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
+                            w, z0 )
+
+       USE control_parameters
+       USE grid_variables
+       USE indices
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i             !< running index x direction
-       INTEGER(iwp) ::  j             !< running index y direction
-       INTEGER(iwp) ::  k             !< running index z direction
-       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
-       INTEGER(iwp) ::  m             !< running index surface elements
-       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
-       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
-
-       REAL(wp) ::  flag              !< flag to mask topography grid points
-       REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
-       REAL(wp) ::  kmxp              !<diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
-       REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
-       REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
-       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
-       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
-       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
-       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
-
-
-
-       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
-       !$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmym, kmyp) &
-       !$ACC PRIVATE(mask_west, mask_east, mask_south, mask_north) &
-       !$ACC PRESENT(wall_flags_total_0, km) &
-       !$ACC PRESENT(u, v, w) &
-       !$ACC PRESENT(ddzu, ddzw, rho_air, drho_air_zw) &
-       !$ACC PRESENT(surf_def_v(0:3)) &
-       !$ACC PRESENT(surf_lsm_v(0:3)) &
-       !$ACC PRESENT(surf_usm_v(0:3)) &
-       !$ACC PRESENT(tend)
+       INTEGER ::  i, j, k
+       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
+                   kmyp_z, wsus, wsvs
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+1), &
+                   km_damp_y(nys-1:nyn+1)
+       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
+
+
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb+1, nzt-1
-!
-!--             Predetermine flag to mask topography and wall-bounded grid points.
-                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) )
-                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
-                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
-                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
-                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
+             DO  k = nzb_w_outer(j,i)+1, nzt-1
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +                                   &
-                                   km(k+1,j,i) + km(k+1,j,i+1) )
-                kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +                                   &
-                                   km(k+1,j,i) + km(k+1,j,i-1) )
-                kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +                                   &
-                                   km(k,j+1,i) + km(k+1,j+1,i) )
-                kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +                                   &
-                                   km(k,j-1,i) + km(k+1,j-1,i) )
-
-                tend(k,j,i) = tend(k,j,i)                                                          &
-                              + ( mask_east *  kmxp * (                                            &
-                                          ( w(k,j,i+1)   - w(k,j,i)   ) * ddx                      &
-                                        + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)                &
-                                                      )                                            &
-                                - mask_west * kmxm *  (                                            &
-                                          ( w(k,j,i)     - w(k,j,i-1) ) * ddx                      &
-                                        + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)                &
-                                                      )                                            &
-                                ) * ddx                                 * flag                     &
-                              + ( mask_north * kmyp * (                                            &
-                                          ( w(k,j+1,i)   - w(k,j,i)   ) * ddy                      &
-                                        + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)                &
-                                                      )                                            &
-                                - mask_south * kmym * (                                            &
-                                          ( w(k,j,i)     - w(k,j-1,i) ) * ddy                      &
-                                        + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)                &
-                                                      )                                            &
-                                ) * ddy                                 * flag                     &
-                              + 2.0_wp * (                                                         &
-                                km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1)              &
-                                            * rho_air(k+1)                                         &
-                              - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)                &
-                                            * rho_air(k)                                           &
-                                         ) * ddzu(k+1) * drho_air_zw(k) * flag
+                kmxp_x = 0.25 * &
+                         ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+                kmxm_x = 0.25 * &
+                         ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+                kmxp_z = kmxp_x
+                kmxm_z = kmxm_x
+                kmyp_y = 0.25 * &
+                         ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+                kmym_y = 0.25 * &
+                         ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                kmyp_z = kmyp_y
+                kmym_z = kmym_y
+!
+!--             Increase diffusion at the outflow boundary in case of
+!--             non-cyclic lateral boundaries. Damping is only needed for
+!--             velocity components parallel to the outflow boundary in
+!--             the direction normal to the outflow boundary.
+                IF ( bc_lr /= 'cyclic' )  THEN
+                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
+                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
+                ENDIF
+                IF ( bc_ns /= 'cyclic' )  THEN
+                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
+                   kmym_y = MAX( kmym_y, km_damp_y(j) )
+                ENDIF
+
+                tend(k,j,i) = tend(k,j,i)                                      &
+                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
+                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
+                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
+                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      &   ) * ddx                                              &
+                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
+                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
+                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
+                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      &   ) * ddy                                              &
+                      & + 2.0 * (                                              &
+                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
+                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
+                      &         ) * ddzu(k+1)
              ENDDO
 
 !
-!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces.
-!--          Note, in the the flat case, loops won't be entered as start_index > end_index.
-!--          Furtermore, note, no vertical natural surfaces so far.
-!--          Default-type surfaces
-             DO  l = 0, 1
-                surf_s = surf_def_v(l)%start_index(j,i)
-                surf_e = surf_def_v(l)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k           = surf_def_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
+!--          Wall functions at all vertical walls, where necessary
+             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
+                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
+                   IF ( wall_w_x(j,i) /= 0.0 )  THEN
+                      wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
+                      wsus = -wsus * ABS( wsus )
+                   ELSE
+                      wsus = 0.0
+                   ENDIF
+                   IF ( wall_w_y(j,i) /= 0.0 )  THEN
+                      wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
+                      wsvs = -wsvs * ABS( wsvs )
+                   ELSE
+                      wsvs = 0.0
+                   ENDIF
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmxp_x = 0.25 * &
+                            ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+                   kmxm_x = 0.25 * &
+                            ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+                   kmxp_z = kmxp_x
+                   kmxm_z = kmxm_x
+                   kmyp_y = 0.25 * &
+                            ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+                   kmym_y = 0.25 * &
+                            ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                   kmyp_z = kmyp_y
+                   kmym_z = kmym_y
+!
+!--                Increase diffusion at the outflow boundary in case of
+!--                non-cyclic lateral boundaries. Damping is only needed for
+!--                velocity components parallel to the outflow boundary in
+!--                the direction normal to the outflow boundary.
+                   IF ( bc_lr /= 'cyclic' )  THEN
+                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
+                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
+                   ENDIF
+                   IF ( bc_ns /= 'cyclic' )  THEN
+                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
+                      kmym_y = MAX( kmym_y, km_damp_y(j) )
+                   ENDIF
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                 + (   fwxp(j,i) * (                           &
+                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
+                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
+                                                   )                           &
+                                     - fwxm(j,i) * (                           &
+                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
+                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
+                                                   )                           &
+                                     + wall_w_x(j,i) * wsus                    &
+                                   ) * ddx                                     &
+                                 + (   fwyp(j,i) * (                           &
+                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
+                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
+                                                   )                           &
+                                     - fwym(j,i) * (                           &
+                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
+                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
+                                                   )                           &
+                                     + wall_w_y(j,i) * wsvs                    &
+                                   ) * ddy                                     &
+                                 + 2.0 * (                                     &
+                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
+                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
+                                         ) * ddzu(k+1)
                 ENDDO
-             ENDDO
-!
-!--          Natural-type surfaces
-             DO  l = 0, 1
-                surf_s = surf_lsm_v(l)%start_index(j,i)
-                surf_e = surf_lsm_v(l)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k           = surf_lsm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
-                ENDDO
-             ENDDO
-!
-!--          Urban-type surfaces
-             DO  l = 0, 1
-                surf_s = surf_usm_v(l)%start_index(j,i)
-                surf_e = surf_usm_v(l)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k           = surf_usm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
-                ENDDO
-             ENDDO
-!
-!--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surface.
-!--          Default-type surfaces
-             DO  l = 2, 3
-                surf_s = surf_def_v(l)%start_index(j,i)
-                surf_e = surf_def_v(l)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k           = surf_def_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
-                ENDDO
-             ENDDO
-!
-!--          Natural-type surfaces
-             DO  l = 2, 3
-                surf_s = surf_lsm_v(l)%start_index(j,i)
-                surf_e = surf_lsm_v(l)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k           = surf_lsm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
-                ENDDO
-             ENDDO
-!
-!--          Urban-type surfaces
-             DO  l = 2, 3
-                surf_s = surf_usm_v(l)%start_index(j,i)
-                surf_e = surf_usm_v(l)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k           = surf_usm_v(l)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
-                ENDDO
-             ENDDO
+             ENDIF
 
           ENDDO
@@ -233,153 +218,135 @@
 
 
-!--------------------------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Call for grid point i,j
-!--------------------------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_w_ij( i, j )
-
-       USE arrays_3d,                                                                              &
-           ONLY :  ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
-
-       USE grid_variables,                                                                         &
-           ONLY :  ddx, ddy
-
-       USE indices,                                                                                &
-           ONLY :  nzb, nzt, wall_flags_total_0
-
-       USE kinds
-
-       USE surface_mod,                                                                            &
-           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
+                               tend, u, v, w, z0 )
+
+       USE control_parameters
+       USE grid_variables
+       USE indices
 
        IMPLICIT NONE
 
-
-       INTEGER(iwp) ::  i             !< running index x direction
-       INTEGER(iwp) ::  j             !< running index y direction
-       INTEGER(iwp) ::  k             !< running index z direction
-       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
-       INTEGER(iwp) ::  m             !< running index surface elements
-       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
-       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
-
-       REAL(wp) ::  flag              !< flag to mask topography grid points
-       REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
-       REAL(wp) ::  kmxp              !< diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
-       REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
-       REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
-       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
-       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
-       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
-       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
-
-
-       DO  k = nzb+1, nzt-1
-!
-!--       Predetermine flag to mask topography and wall-bounded grid points.
-          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) )
-          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
-          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
-          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
-          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
+       INTEGER ::  i, j, k
+       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
+                   kmyp_z, wsus, wsvs
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+1), &
+                   km_damp_y(nys-1:nyn+1)
+       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
+
+
+       DO  k = nzb_w_outer(j,i)+1, nzt-1
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-          kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-          kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-
-          tend(k,j,i) = tend(k,j,i)                                                                &
-                        + ( mask_east *  kmxp * (                                                  &
-                                    ( w(k,j,i+1)   - w(k,j,i)   ) * ddx                            &
-                                  + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)                      &
-                                                )                                                  &
-                          - mask_west * kmxm *  (                                                  &
-                                    ( w(k,j,i)     - w(k,j,i-1) ) * ddx                            &
-                                  + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)                      &
-                                                )                                                  &
-                          ) * ddx                                 * flag                           &
-                        + ( mask_north * kmyp * (                                                  &
-                                    ( w(k,j+1,i)   - w(k,j,i)   ) * ddy                            &
-                                  + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)                      &
-                                                )                                                  &
-                          - mask_south * kmym * (                                                  &
-                                    ( w(k,j,i)     - w(k,j-1,i) ) * ddy                            &
-                                  + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)                      &
-                                                )                                                  &
-                          ) * ddy                                 * flag                           &
-                        + 2.0_wp * (                                                               &
-                          km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)                      &
-                                      * rho_air(k+1)                                               &
-                        - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)                      &
-                                      * rho_air(k)                                                 &
-                                   ) * ddzu(k+1) * drho_air_zw(k) * flag
+          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+          kmxp_z = kmxp_x
+          kmxm_z = kmxm_x
+          kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+          kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+          kmyp_z = kmyp_y
+          kmym_z = kmym_y
+!
+!--       Increase diffusion at the outflow boundary in case of non-cyclic
+!--       lateral boundaries. Damping is only needed for velocity components
+!--       parallel to the outflow boundary in the direction normal to the
+!--       outflow boundary.
+          IF ( bc_lr /= 'cyclic' )  THEN
+             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
+             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
+          ENDIF
+          IF ( bc_ns /= 'cyclic' )  THEN
+             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
+             kmym_y = MAX( kmym_y, km_damp_y(j) )
+          ENDIF
+
+          tend(k,j,i) = tend(k,j,i)                                            &
+                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
+                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
+                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
+                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      &   ) * ddx                                              &
+                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
+                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
+                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
+                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      &   ) * ddy                                              &
+                      & + 2.0 * (                                              &
+                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
+                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
+                      &         ) * ddzu(k+1)
        ENDDO
-!
-!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces. Note, in
-!--    the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
-!--    vertical natural surfaces so far.
-!--    Default-type surfaces
-       DO  l = 0, 1
-          surf_s = surf_def_v(l)%start_index(j,i)
-          surf_e = surf_def_v(l)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k           = surf_def_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
+
+!
+!--    Wall functions at all vertical walls, where necessary
+       IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
+          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
+             IF ( wall_w_x(j,i) /= 0.0 )  THEN
+                wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
+                wsus = -wsus * ABS( wsus )
+             ELSE
+                wsus = 0.0
+             ENDIF
+             IF ( wall_w_y(j,i) /= 0.0 )  THEN
+                wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
+                wsvs = -wsvs * ABS( wsvs )
+             ELSE
+                wsvs = 0.0
+             ENDIF
+!
+!--          Interpolate eddy diffusivities on staggered gridpoints
+             kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+             kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+             kmxp_z = kmxp_x
+             kmxm_z = kmxm_x
+             kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+             kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+             kmyp_z = kmyp_y
+             kmym_z = kmym_y
+!
+!--          Increase diffusion at the outflow boundary in case of
+!--          non-cyclic lateral boundaries. Damping is only needed for
+!--          velocity components parallel to the outflow boundary in
+!--          the direction normal to the outflow boundary.
+             IF ( bc_lr /= 'cyclic' )  THEN
+                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
+                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
+             ENDIF
+             IF ( bc_ns /= 'cyclic' )  THEN
+                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
+                kmym_y = MAX( kmym_y, km_damp_y(j) )
+             ENDIF
+
+             tend(k,j,i) = tend(k,j,i)                                         &
+                                 + (   fwxp(j,i) * (                           &
+                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
+                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
+                                                   )                           &
+                                     - fwxm(j,i) * (                           &
+                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
+                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
+                                                   )                           &
+                                     + wall_w_x(j,i) * wsus                    &
+                                   ) * ddx                                     &
+                                 + (   fwyp(j,i) * (                           &
+                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
+                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
+                                                   )                           &
+                                     - fwym(j,i) * (                           &
+                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
+                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
+                                                   )                           &
+                                     + wall_w_y(j,i) * wsvs                    &
+                                   ) * ddy                                     &
+                                 + 2.0 * (                                     &
+                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
+                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
+                                         ) * ddzu(k+1)
           ENDDO
-       ENDDO
-!
-!--    Natural-type surfaces
-       DO  l = 0, 1
-          surf_s = surf_lsm_v(l)%start_index(j,i)
-          surf_e = surf_lsm_v(l)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k           = surf_lsm_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
-          ENDDO
-       ENDDO
-!
-!--    Urban-type surfaces
-       DO  l = 0, 1
-          surf_s = surf_usm_v(l)%start_index(j,i)
-          surf_e = surf_usm_v(l)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k           = surf_usm_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
-          ENDDO
-       ENDDO
-!
-!--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surfaces.
-!--    Default-type surfaces
-       DO  l = 2, 3
-          surf_s = surf_def_v(l)%start_index(j,i)
-          surf_e = surf_def_v(l)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k           = surf_def_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
-          ENDDO
-       ENDDO
-!
-!--    Natural-type surfaces
-       DO  l = 2, 3
-          surf_s = surf_lsm_v(l)%start_index(j,i)
-          surf_e = surf_lsm_v(l)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k           = surf_lsm_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
-          ENDDO
-       ENDDO
-!
-!--    Urban-type surfaces
-       DO  l = 2, 3
-          surf_s = surf_usm_v(l)%start_index(j,i)
-          surf_e = surf_usm_v(l)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k           = surf_usm_v(l)%k(m)
-             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
-          ENDDO
-       ENDDO
-
+       ENDIF
 
     END SUBROUTINE diffusion_w_ij
