Index: palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_w.f90	(revision 2119)
+++ palm/trunk/SOURCE/diffusion_w.f90	(revision 2232)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adjustments to new topography and surface concept
 ! 
 ! Former revisions:
@@ -97,7 +97,4 @@
  
 
-    USE wall_fluxes_mod,                                                       &
-        ONLY :  wall_fluxes
-
     PRIVATE
     PUBLIC diffusion_w
@@ -125,114 +122,161 @@
            
        USE grid_variables,                                                     &     
-           ONLY :  ddx, ddy, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
+           ONLY :  ddx, ddy
            
        USE indices,                                                            &            
-           ONLY :  nxl, nxr, nyn, nys, nzb, nzb_w_inner, nzb_w_outer, nzt
+           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
            
        USE kinds
 
+       USE surface_mod,                                                        &
+           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i     !<
-       INTEGER(iwp) ::  j     !<
-       INTEGER(iwp) ::  k     !<
+       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) ::  kmxm  !<
-       REAL(wp) ::  kmxp  !<
-       REAL(wp) ::  kmym  !<
-       REAL(wp) ::  kmyp  !< 
-
-       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !<
-       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !<
-
-
-!
-!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
-!--    walls, if neccessary
-       IF ( topography /= 'flat' )  THEN
-          CALL wall_fluxes( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, nzb_w_inner,             &
-                            nzb_w_outer, wall_w_x )
-          CALL wall_fluxes( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, nzb_w_inner,             &
-                            nzb_w_outer, wall_w_y )
-       ENDIF
+       REAL(wp) ::  flag              !< flag to mask topography grid points
+       REAL(wp) ::  kmxm              !<
+       REAL(wp) ::  kmxp              !<
+       REAL(wp) ::  kmym              !<
+       REAL(wp) ::  kmyp              !< 
+       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
+       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point 
+       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point 
+       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point 
+
+
 
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb_w_outer(j,i)+1, nzt-1
+             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_0(k,j,i),   3 ) ) 
+                mask_east  = MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(k,j,i+1), 3 ) )
+                mask_west  = MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(k,j,i-1), 3 ) )
+                mask_south = MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(k,j-1,i), 3 ) )
+                mask_north = MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(k,j+1,i), 3 ) )
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmxp = 0.25_wp *                                               &
-                       ( km(k,j,i)+km(k,j,i+1)+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) )
+                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)                                      &
-                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
-                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
-                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
-                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
-                      &   ) * ddx                                              &
-                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
-                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
-                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
-                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
-                      &   ) * ddy                                              &
-                      & + 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)
-             ENDDO
-
-!
-!--          Wall functions at all vertical walls, where necessary
-             IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
-
-                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
-!
-!--                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)                                   &
-                                 + (   fwxp(j,i) * (                           &
-                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
-                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
-                                                   )                           &
-                                     - fwxm(j,i) * (                           &
-                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
-                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
-                                                   )                           &
-                                     + wall_w_x(j,i) * wsus(k,j,i)             &
-                                   ) * ddx                                     &
-                                 + (   fwyp(j,i) * (                           &
-                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
-                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
-                                                   )                           &
-                                     - fwym(j,i) * (                           &
-                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
-                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
-                                                   )                           &
-                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
-                                   ) * ddy                                     &
-                                 + 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)
-                ENDDO
-             ENDIF
+                       + ( mask_east *  kmxp * (                               &
+                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
+                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
+                                               )                               &
+                         - mask_west * kmxm *  (                               &
+                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
+                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
+                                               )                               &
+                         ) * ddx                                 * flag        &
+                       + ( mask_north * kmyp * (                               &
+                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
+                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
+                                               )                               &
+                         - mask_south * kmym * (                               &
+                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
+                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
+                                               )                               &
+                         ) * ddy                                 * flag        &
+                       + 2.0_wp * (                                            &
+                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1) &
+                                     * rho_air(k+1)                            &
+                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
+                                     * rho_air(k)                              &
+                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
+             ENDDO
+
+!
+!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
+!--          surfaces. Note, in the the flat case, loops won't be entered as 
+!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
+!--          so far.           
+!--          Default-type surfaces
+             DO  l = 0, 1
+                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
+                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
 
           ENDDO
@@ -256,27 +300,44 @@
            
        USE grid_variables,                                                     &     
-           ONLY :  ddx, ddy, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
+           ONLY :  ddx, ddy
            
        USE indices,                                                            &            
-           ONLY :  nxl, nxr, nyn, nys, nzb, nzb_w_inner, nzb_w_outer, nzt
+           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
            
        USE kinds
 
+       USE surface_mod,                                                        &
+           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i     !<
-       INTEGER(iwp) ::  j     !<
-       INTEGER(iwp) ::  k     !<
+
+       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) ::  kmxm  !<
-       REAL(wp) ::  kmxp  !<
-       REAL(wp) ::  kmym  !<
-       REAL(wp) ::  kmyp  !< 
-
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs
-
-
-       DO  k = nzb_w_outer(j,i)+1, nzt-1
+       REAL(wp) ::  flag              !< flag to mask topography grid points
+       REAL(wp) ::  kmxm              !<
+       REAL(wp) ::  kmxp              !<
+       REAL(wp) ::  kmym              !<
+       REAL(wp) ::  kmyp              !< 
+       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
+       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point 
+       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point 
+       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point 
+
+
+       DO  k = nzb+1, nzt-1
+!
+!--       Predetermine flag to mask topography and wall-bounded grid points. 
+          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   3 ) ) 
+          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 3 ) )
+          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 3 ) )
+          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 3 ) )
+          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 3 ) )
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
@@ -287,79 +348,102 @@
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
-                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
-                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
-                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
-                      &   ) * ddx                                              &
-                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
-                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
-                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
-                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
-                      &   ) * ddy                                              &
-                      & + 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)
-       ENDDO
-
-!
-!--    Wall functions at all vertical walls, where necessary
-       IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
-
-!
-!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
-          IF ( wall_w_x(j,i) /= 0.0_wp )  THEN
-             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
-                               wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
-          ELSE
-             wsus = 0.0_wp
-          ENDIF
-
-          IF ( wall_w_y(j,i) /= 0.0_wp )  THEN
-             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
-                               wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
-          ELSE
-             wsvs = 0.0_wp
-          ENDIF
-
-          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
-!
-!--          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)                                         &
-                                 + (   fwxp(j,i) * (                           &
-                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
-                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
-                                                   )                           &
-                                     - fwxm(j,i) * (                           &
-                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
-                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
-                                                   )                           &
-                                     + wall_w_x(j,i) * wsus(k)                 &
-                                   ) * ddx                                     &
-                                 + (   fwyp(j,i) * (                           &
-                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
-                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
-                                                   )                           &
-                                     - fwym(j,i) * (                           &
-                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
-                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
-                                                   )                           &
-                                     + wall_w_y(j,i) * wsvs(k)                 &
-                                   ) * ddy                                     &
-                                 + 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)
-          ENDDO
-       ENDIF
+                       + ( mask_east *  kmxp * (                               &
+                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
+                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
+                                               )                               &
+                         - mask_west * kmxm *  (                               &
+                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
+                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
+                                               )                               &
+                         ) * ddx                                 * flag        &
+                       + ( mask_north * kmyp * (                               &
+                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
+                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
+                                               )                               &
+                         - mask_south * kmym * (                               &
+                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
+                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
+                                               )                               &
+                         ) * ddy                                 * flag        &
+                       + 2.0_wp * (                                            &
+                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)   &
+                                     * rho_air(k+1)                            &
+                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
+                                     * rho_air(k)                              &
+                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
+       ENDDO
+!
+!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
+!--    surfaces. Note, in the the flat case, loops won't be entered as 
+!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
+!--    so far.           
+!--    Default-type surfaces
+       DO  l = 0, 1
+          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
+          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
+
 
     END SUBROUTINE diffusion_w_ij
