Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 2127)
+++ palm/trunk/SOURCE/production_e.f90	(revision 2232)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adjustments to new surface concept
 !
 ! Former revisions:
@@ -106,8 +106,4 @@
 !------------------------------------------------------------------------------!
  MODULE production_e_mod
- 
-
-    USE wall_fluxes_mod,                                                       &
-        ONLY:  wall_fluxes_e
 
     USE kinds
@@ -115,9 +111,4 @@
     PRIVATE
     PUBLIC production_e, production_e_init
-
-    LOGICAL, SAVE ::  first_call = .TRUE.  !<
-
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  u_0  !<
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  v_0  !<
 
     INTERFACE production_e
@@ -141,6 +132,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, qsws, qswst, shf,      &
-                  tend, tswst, u, v, vpt, w
+           ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w
 
        USE cloud_parameters,                                                   &
@@ -154,88 +144,79 @@
 
        USE grid_variables,                                                     &
-           ONLY:  ddx, dx, ddy, dy, wall_e_x, wall_e_y
+           ONLY:  ddx, dx, ddy, dy
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner,                   &
-                   nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
+
+       USE surface_mod,                                                        &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
+                   surf_usm_v
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i           !<
-       INTEGER(iwp) ::  j           !<
-       INTEGER(iwp) ::  k           !<
-
-       REAL(wp)     ::  def         !<
-       REAL(wp)     ::  dudx        !<
-       REAL(wp)     ::  dudy        !<
-       REAL(wp)     ::  dudz        !<
-       REAL(wp)     ::  dvdx        !<
-       REAL(wp)     ::  dvdy        !<
-       REAL(wp)     ::  dvdz        !<
-       REAL(wp)     ::  dwdx        !<
-       REAL(wp)     ::  dwdy        !<
-       REAL(wp)     ::  dwdz        !<
+       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 for different surface type orientation
+       INTEGER(iwp) ::  m       !< running index surface elements
+       INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
+       INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
+
+       REAL(wp)     ::  def         !< 
+       REAL(wp)     ::  flag        !< flag to mask topography
        REAL(wp)     ::  k1          !<
        REAL(wp)     ::  k2          !<
-       REAL(wp)     ::  km_neutral  !<
+       REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
        REAL(wp)     ::  theta       !<
        REAL(wp)     ::  temp        !<
-
-!       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  usvs  !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus  !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs  !<
-
-!
-!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
-!--    vertical walls, if neccessary
-!--    So far, results are slightly different from the ij-Version.
-!--    Therefore, ij-Version is called further below within the ij-loops.
-!       IF ( topography /= 'flat' )  THEN
-!          CALL wall_fluxes_e( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, wall_e_y )
-!          CALL wall_fluxes_e( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, wall_e_y )
-!          CALL wall_fluxes_e( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, wall_e_x )
-!          CALL wall_fluxes_e( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, wall_e_x )
-!       ENDIF
-
+       REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation 
+       REAL(wp)     ::  usvs        !< momentum flux u"v"
+       REAL(wp)     ::  vsus        !< momentum flux v"u"
+       REAL(wp)     ::  wsus        !< momentum flux w"u"
+       REAL(wp)     ::  wsvs        !< momentum flux w"v"
+
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dudx  !< Gradient of u-component in x-direction
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dudy  !< Gradient of u-component in y-direction
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dudz  !< Gradient of u-component in z-direction
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dvdx  !< Gradient of v-component in x-direction
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dvdy  !< Gradient of v-component in y-direction
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dvdz  !< Gradient of v-component in z-direction
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dwdx  !< Gradient of w-component in x-direction
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dwdy  !< Gradient of w-component in y-direction
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dwdz  !< Gradient of w-component in z-direction
 
        DO  i = nxl, nxr
 
-!
-!--       Calculate TKE production by shear
-          DO  j = nys, nyn
-             DO  k = nzb_diff_s_outer(j,i), nzt
-
-                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0_wp )  def = 0.0_wp
-
-                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
+          IF ( constant_flux_layer )  THEN
+
+!
+!--          Calculate TKE production by shear. Calculate gradients at all grid 
+!--          points first, gradients at surface-bounded grid points will be 
+!--          overwritten further below. 
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+
+                   dudx(k,j) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -         &
+                                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dudz(k,j) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -         &
+                                           u(k-1,j,i) - u(k-1,j,i+1) ) *       &
+                                                            dd2zu(k)
+   
+                   dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -         &
+                                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy(k,j) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz(k,j) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -         &
+                                           v(k-1,j,i) - v(k-1,j+1,i) ) *       &
+                                                            dd2zu(k)
+
+                   dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -         &
+                                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -         &
+                                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dwdz(k,j) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                ENDDO
              ENDDO
-          ENDDO
-
-          IF ( constant_flux_layer )  THEN
 
 !
@@ -244,282 +225,251 @@
 !--          'bottom and wall: use u_0,v_0 and wall functions'
              DO  j = nys, nyn
-
-                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
-                THEN
-
-                   k = nzb_diff_s_inner(j,i) - 1
-                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                   dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                     u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
-                   dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                     v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
-
-                   IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
-!
-!--                   Inconsistency removed: as the thermal stratification is 
-!--                   not taken into account for the evaluation of the wall 
-!--                   fluxes at vertical walls, the eddy viscosity km must not 
-!--                   be used for the evaluation of the velocity gradients dudy 
-!--                   and dwdy
-!--                   Note: The validity of the new method has not yet been 
-!--                         shown, as so far no suitable data for a validation 
-!--                         has been available
-                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
-                                          usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
-                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
-                                          wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
-                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
-                                   0.5_wp * dy 
-                      IF ( km_neutral > 0.0_wp )  THEN
-                         dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
-                         dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
-                      ELSE
-                         dudy = 0.0_wp
-                         dwdy = 0.0_wp
-                      ENDIF
-                   ELSE
-                      dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                      dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   ENDIF
-
-                   IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
-!
-!--                   Inconsistency removed: as the thermal stratification is 
-!--                   not taken into account for the evaluation of the wall 
-!--                   fluxes at vertical walls, the eddy viscosity km must not 
-!--                   be used for the evaluation of the velocity gradients dvdx 
-!--                   and dwdx
-!--                   Note: The validity of the new method has not yet been 
-!--                         shown, as so far no suitable data for a validation 
-!--                         has been available
-                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
-                                          vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
-                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
-                                          wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
-                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * &
-                                   0.5_wp * dx
-                      IF ( km_neutral > 0.0_wp )  THEN
-                         dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
-                         dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
-                      ELSE
-                         dvdx = 0.0_wp
-                         dwdx = 0.0_wp
-                      ENDIF
-                   ELSE
-                      dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                      dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   ENDIF
-
-                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+!
+!--             Compute gradients at north- and south-facing surfaces.
+!--             First, for default surfaces, then for urban surfaces. 
+!--             Note, so far no natural vertical surfaces implemented
+                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)
+                      usvs        = surf_def_v(l)%mom_flux_tke(0,m)
+                      wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
+      
+                      km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp      &
+                                      * 0.5_wp * dy
+!
+!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
+                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
+                                        BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+                      dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
+                      dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
+                   ENDDO
+!
+!--                Natural surfaces
+                   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)
+                      usvs        = surf_lsm_v(l)%mom_flux_tke(0,m)
+                      wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
+      
+                      km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp      &
+                                      * 0.5_wp * dy
+!
+!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
+                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
+                                        BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+                      dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
+                      dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
+                   ENDDO 
+!
+!--                Urban surfaces
+                   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)
+                      usvs        = surf_usm_v(l)%mom_flux_tke(0,m)
+                      wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
+      
+                      km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp      &
+                                      * 0.5_wp * dy
+!
+!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
+                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
+                                        BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+                      dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
+                      dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
+                   ENDDO  
+                ENDDO
+!
+!--             Compute gradients at east- and west-facing walls
+                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)
+                      vsus  = surf_def_v(l)%mom_flux_tke(0,m)
+                      wsus  = surf_def_v(l)%mom_flux_tke(1,m)
+
+                      km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp      &
+                                         * 0.5_wp * dx
+!
+!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
+                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
+                                        BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+                      dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
+                      dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
+                   ENDDO 
+!
+!--                Natural surfaces                   
+                   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)
+                      vsus  = surf_lsm_v(l)%mom_flux_tke(0,m)
+                      wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
+
+                      km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp      &
+                                         * 0.5_wp * dx
+!
+!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
+                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
+                                        BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+                      dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
+                      dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
+                   ENDDO   
+!
+!--                Urban surfaces                   
+                   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)
+                      vsus  = surf_usm_v(l)%mom_flux_tke(0,m)
+                      wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
+
+                      km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp      &
+                                         * 0.5_wp * dx
+!
+!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
+                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
+                                        BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+                      dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
+                      dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
+                   ENDDO 
+                ENDDO
+!
+!--             Compute gradients at upward-facing surfaces
+                surf_s = surf_def_h(0)%start_index(j,i)
+                surf_e = surf_def_h(0)%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_def_h(0)%k(m)
+!
+!--                Please note, actually, an interpolation of u_0 and v_0 
+!--                onto the grid center would be required. However, this 
+!--                would require several data transfers between 2D-grid and
+!--                wall type. The effect of this missing interpolation is 
+!--                negligible. (See also production_e_init).
+                   dudz(k,j) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)   
+                   dvdz(k,j) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
+            
+                ENDDO
+!
+!--             Natural surfaces
+                surf_s = surf_lsm_h%start_index(j,i)
+                surf_e = surf_lsm_h%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_lsm_h%k(m)
+!
+!--                Please note, actually, an interpolation of u_0 and v_0 
+!--                onto the grid center would be required. However, this 
+!--                would require several data transfers between 2D-grid and
+!--                wall type. The effect of this missing interpolation is 
+!--                negligible. (See also production_e_init).
+                   dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
+                   dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
+            
+                ENDDO
+!
+!--             Urban surfaces
+                surf_s = surf_usm_h%start_index(j,i)
+                surf_e = surf_usm_h%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_usm_h%k(m)
+!
+!--                Please note, actually, an interpolation of u_0 and v_0 
+!--                onto the grid center would be required. However, this 
+!--                would require several data transfers between 2D-grid and
+!--                wall type. The effect of this missing interpolation is 
+!--                negligible. (See also production_e_init).
+                   dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
+                   dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
+            
+                ENDDO
+!
+!--             Compute gradients at downward-facing walls, only for 
+!--             non-natural default surfaces
+                surf_s = surf_def_h(1)%start_index(j,i)
+                surf_e = surf_def_h(1)%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_def_h(1)%k(m)
+!
+!--                Please note, actually, an interpolation of u_0 and v_0 
+!--                onto the grid center would be required. However, this 
+!--                would require several data transfers between 2D-grid and
+!--                wall type. The effect of this missing interpolation is 
+!--                negligible. (See also production_e_init).
+                   dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
+                   dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
+
+                ENDDO
+
+             ENDDO
+
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+
+                   def = 2.0_wp * ( dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 ) + &
+                                    dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 +   &
+                                    dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 +   &
+                         2.0_wp * ( dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j)  +   &
+                                    dwdy(k,j)*dvdz(k,j) )
 
                    IF ( def < 0.0_wp )  def = 0.0_wp
 
-                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-
-!
-!--                (3) - will be executed only, if there is at least one level
-!--                between (2) and (4), i.e. the topography must have a
-!--                minimum height of 2 dz. Wall fluxes for this case have
-!--                already been calculated for (2).
-!--                'wall only: use wall functions'
-
-                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
-
-                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                      dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-                      dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                      dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
-
-                      IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
-!
-!--                      Inconsistency removed: as the thermal stratification 
-!--                      is not taken into account for the evaluation of the 
-!--                      wall fluxes at vertical walls, the eddy viscosity km 
-!--                      must not be used for the evaluation of the velocity 
-!--                      gradients dudy and dwdy
-!--                      Note: The validity of the new method has not yet 
-!--                            been shown, as so far no suitable data for a 
-!--                            validation has been available
-                         km_neutral = kappa * ( usvs(k)**2 + & 
-                                                wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
-                         IF ( km_neutral > 0.0_wp )  THEN
-                            dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
-                            dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
-                         ELSE
-                            dudy = 0.0_wp
-                            dwdy = 0.0_wp
-                         ENDIF
-                      ELSE
-                         dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                   flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
+
+                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
+
+                ENDDO
+             ENDDO
+
+          ELSEIF ( use_surface_fluxes )  THEN
+
+             DO  j = nys, nyn
+!
+!--             Calculate TKE production by shear. Here, no additional
+!--             wall-bounded code is considered. 
+!--             Why?
+                DO  k = nzb+1, nzt
+
+                   dudx(k,j)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy(k,j)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -        &
                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                         dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                   dudz(k,j)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -        &
+                                            u(k-1,j,i) - u(k-1,j,i+1) ) *      &
+                                                                        dd2zu(k)
+
+                   dvdx(k,j)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -        &
+                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy(k,j)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz(k,j)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -        &
+                                            v(k-1,j,i) - v(k-1,j+1,i) ) *      &
+                                                                        dd2zu(k)
+
+                   dwdx(k,j)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -        &
+                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy(k,j)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -        &
                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                      ENDIF
-
-                      IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
-!
-!--                      Inconsistency removed: as the thermal stratification 
-!--                      is not taken into account for the evaluation of the 
-!--                      wall fluxes at vertical walls, the eddy viscosity km 
-!--                      must not be used for the evaluation of the velocity 
-!--                      gradients dvdx and dwdx
-!--                      Note: The validity of the new method has not yet 
-!--                            been shown, as so far no suitable data for a 
-!--                            validation has been available
-                         km_neutral = kappa * ( vsus(k)**2 + & 
-                                                wsus(k)**2 )**0.25_wp * 0.5_wp * dx
-                         IF ( km_neutral > 0.0_wp )  THEN
-                            dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
-                            dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
-                         ELSE
-                            dvdx = 0.0_wp
-                            dwdx = 0.0_wp
-                         ENDIF
-                      ELSE
-                         dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                         dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                      ENDIF
-
-                      def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                           dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
-                           dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                      IF ( def < 0.0_wp )  def = 0.0_wp
-
-                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-                   ENDDO
-
-                ENDIF
-
-             ENDDO
-
-!
-!--          (4) - will allways be executed.
-!--          'special case: free atmosphere' (as for case (0))
-             DO  j = nys, nyn
-
-                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
-                THEN
-
-                   k = nzb_diff_s_outer(j,i)-1
-
-                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+                   dwdz(k,j)  =           ( w(k,j,i)   - w(k-1,j,i)   ) *      &
+                                                                        ddzw(k)
+      
+                   def = 2.0_wp * (                                            &
+                                dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2     &
+                                  ) +                                          &
+                                dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 +   &
+                                dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 +   &
+                         2.0_wp * (                                            &
+                                dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j)  +   &
+                                dwdy(k,j)*dvdz(k,j)                            &
+                                  )
 
                    IF ( def < 0.0_wp )  def = 0.0_wp
 
-                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-                ENDIF
-
-             ENDDO
-
-!
-!--          Position without adjacent wall
-!--          (1) - will allways be executed. 
-!--          'bottom only: use u_0,v_0'
-             DO  j = nys, nyn
-
-                IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
-                THEN
-
-                   k = nzb_diff_s_inner(j,i)-1
-
-                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                       u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-
-                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                       v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0_wp )  def = 0.0_wp
-
-                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-                ENDIF
-
-             ENDDO
-
-          ELSEIF ( use_surface_fluxes )  THEN
-
-             DO  j = nys, nyn
-
-                k = nzb_diff_s_outer(j,i)-1
-
-                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                   u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0_wp )  def = 0.0_wp
-
-                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
+                   flag  = MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,i), 29 ) ) 
+                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
+      
+                ENDDO
              ENDDO
 
@@ -539,9 +489,12 @@
 !--                   in the bottom and top surface layer
                       DO  j = nys, nyn
-                         DO  k = nzb_s_inner(j,i)+1, nzt
+                         DO  k = nzb+1, nzt
                             tend(k,j,i) = tend(k,j,i) +                        &
                                           kh(k,j,i) * g / rho_reference *      &
                                           ( prho(k+1,j,i) - prho(k-1,j,i) ) *  &
-                                          dd2zu(k)
+                                          dd2zu(k) *                           &
+                                    MERGE( 1.0_wp, 0.0_wp,                     &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                          )
                          ENDDO
                       ENDDO
@@ -550,21 +503,60 @@
 
                       DO  j = nys, nyn
-                         DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                            tend(k,j,i) = tend(k,j,i) -                   &
-                                          kh(k,j,i) * g / pt_reference *  &
-                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) * &
-                                          dd2zu(k)
+                         DO  k = nzb+1, nzt
+!
+!--                         Flag 9 is used to mask top fluxes, flag 30 to mask
+!--                         surface fluxes
+                            tend(k,j,i) = tend(k,j,i) -                        &
+                                          kh(k,j,i) * g / pt_reference  *      &
+                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) *      &
+                                          dd2zu(k)                      *      &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           )                            *      &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )  
                          ENDDO
 
                          IF ( use_surface_fluxes )  THEN
-                            k = nzb_diff_s_inner(j,i)-1
-                            tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
-                                                        shf(j,i)
+!
+!--                         Default surfaces, up- and downward-facing
+                            DO  l = 0, 1
+                               surf_s = surf_def_h(l)%start_index(j,i)
+                               surf_e = surf_def_h(l)%end_index(j,i)
+                               DO  m = surf_s, surf_e
+                                  k = surf_def_h(l)%k(m)
+                                  tend(k,j,i) = tend(k,j,i) + g / pt_reference &
+                                                         * surf_def_h(l)%shf(m)   
+                               ENDDO   
+                            ENDDO      
+!
+!--                         Natural surfaces
+                            surf_s = surf_lsm_h%start_index(j,i)
+                            surf_e = surf_lsm_h%end_index(j,i)
+                            DO  m = surf_s, surf_e
+                               k = surf_lsm_h%k(m)
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * surf_lsm_h%shf(m)   
+                            ENDDO
+!
+!--                         Urban surfaces
+                            surf_s = surf_usm_h%start_index(j,i)
+                            surf_e = surf_usm_h%end_index(j,i)
+                            DO  m = surf_s, surf_e
+                               k = surf_usm_h%k(m)
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * surf_usm_h%shf(m)   
+                            ENDDO                          
                          ENDIF
 
                          IF ( use_top_fluxes )  THEN
-                            k = nzt
-                            tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
-                                                        tswst(j,i)
+                            surf_s = surf_def_h(2)%start_index(j,i)
+                            surf_e = surf_def_h(2)%end_index(j,i)
+                            DO  m = surf_s, surf_e
+                               k = surf_def_h(2)%k(m)
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
+                                                           surf_def_h(2)%shf(m) 
+                            ENDDO
                          ENDIF
                       ENDDO
@@ -579,9 +571,12 @@
 !--                   in the bottom and top surface layer
                       DO  j = nys, nyn
-                         DO  k = nzb_s_inner(j,i)+1, nzt
+                         DO  k = nzb+1, nzt
                             tend(k,j,i) = tend(k,j,i) +                        &
                                           kh(k,j,i) * g / prho(k,j,i) *        &
                                           ( prho(k+1,j,i) - prho(k-1,j,i) ) *  &
-                                          dd2zu(k)
+                                          dd2zu(k)                           * &
+                                     MERGE( 1.0_wp, 0.0_wp,                    &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                          )
                          ENDDO
                       ENDDO
@@ -590,21 +585,60 @@
 
                       DO  j = nys, nyn
-                         DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                            tend(k,j,i) = tend(k,j,i) -                   &
-                                          kh(k,j,i) * g / pt(k,j,i) *     &
-                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) * &
-                                          dd2zu(k)
+                         DO  k = nzb+1, nzt
+!
+!--                         Flag 9 is used to mask top fluxes, flag 30 to mask
+!--                         surface fluxes
+                            tend(k,j,i) = tend(k,j,i) -                        &
+                                          kh(k,j,i) * g / pt(k,j,i) *          &
+                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) *      &
+                                          dd2zu(k)                      *      &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           )                            *      &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           ) 
                          ENDDO
 
                          IF ( use_surface_fluxes )  THEN
-                            k = nzb_diff_s_inner(j,i)-1
-                            tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
-                                                        shf(j,i)
+!
+!--                         Default surfaces, up- and downwrd-facing
+                            DO  l = 0, 1
+                               surf_s = surf_def_h(l)%start_index(j,i)
+                               surf_e = surf_def_h(l)%end_index(j,i)
+                               DO  m = surf_s, surf_e
+                                  k = surf_def_h(l)%k(m)
+                                  tend(k,j,i) = tend(k,j,i) + g / pt_reference &
+                                                         * surf_def_h(l)%shf(m)   
+                               ENDDO  
+                            ENDDO
+!
+!--                         Natural surfaces
+                            surf_s = surf_lsm_h%start_index(j,i)
+                            surf_e = surf_lsm_h%end_index(j,i)
+                            DO  m = surf_s, surf_e
+                               k = surf_lsm_h%k(m)
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * surf_lsm_h%shf(m)   
+                            ENDDO  
+!
+!--                         Urban surfaces
+                            surf_s = surf_usm_h%start_index(j,i)
+                            surf_e = surf_usm_h%end_index(j,i)
+                            DO  m = surf_s, surf_e
+                               k = surf_usm_h%k(m)
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * surf_usm_h%shf(m)   
+                            ENDDO  
                          ENDIF
 
                          IF ( use_top_fluxes )  THEN
-                            k = nzt
-                            tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
-                                                        tswst(j,i)
+                            surf_s = surf_def_h(2)%start_index(j,i)
+                            surf_e = surf_def_h(2)%end_index(j,i)
+                            DO  m = surf_s, surf_e
+                               k = surf_def_h(2)%k(m)
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
+                                                           surf_def_h(2)%shf(m) 
+                            ENDDO
                          ENDIF
                       ENDDO
@@ -618,6 +652,8 @@
                 DO  j = nys, nyn
 
-                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
-
+                   DO  k = nzb+1, nzt
+!
+!--                   Flag 9 is used to mask top fluxes, flag 30 to mask
+!--                   surface fluxes
                       IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
@@ -627,48 +663,11 @@
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
                                            k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
-                                         ) * dd2zu(k)
-                      ELSE IF ( cloud_physics )  THEN
-                         IF ( ql(k,j,i) == 0.0_wp )  THEN
-                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                            k2 = 0.61_wp * pt(k,j,i)
-                         ELSE
-                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                            temp  = theta * t_d_pt(k)
-                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
-                                          ( q(k,j,i) - ql(k,j,i) ) *          &
-                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
-                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
-                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
-                         ENDIF
-                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
-                                         g / vpt(k,j,i) *                      &
-                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
-                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
-                                         ) * dd2zu(k)
-                      ELSE IF ( cloud_droplets )  THEN
-                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
-                         k2 = 0.61_wp * pt(k,j,i)
-                         tend(k,j,i) = tend(k,j,i) -                          & 
-                                       kh(k,j,i) * g / vpt(k,j,i) *           &
-                                       ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +  &
-                                         k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -  &
-                                         pt(k,j,i) * ( ql(k+1,j,i) -          &
-                                         ql(k-1,j,i) ) ) * dd2zu(k)
-                      ENDIF
-
-                   ENDDO
-
-                ENDDO
-
-                IF ( use_surface_fluxes )  THEN
-
-                   DO  j = nys, nyn
-
-                      k = nzb_diff_s_inner(j,i)-1
-
-                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                         k2 = 0.61_wp * pt(k,j,i)
+                                         ) * dd2zu(k) *                        &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           )          *                        &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )
                       ELSE IF ( cloud_physics )  THEN
                          IF ( ql(k,j,i) == 0.0_wp )  THEN
@@ -685,13 +684,143 @@
                             k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
+                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
+                                         g / vpt(k,j,i) *                      &
+                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
+                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
+                                         ) * dd2zu(k) *                        &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           )          *                        &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )
                       ELSE IF ( cloud_droplets )  THEN
                          k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
                          k2 = 0.61_wp * pt(k,j,i)
+                         tend(k,j,i) = tend(k,j,i) -                           & 
+                                       kh(k,j,i) * g / vpt(k,j,i) *            &
+                                       ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +   &
+                                         k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -   &
+                                         pt(k,j,i) * ( ql(k+1,j,i) -           &
+                                         ql(k-1,j,i) ) ) * dd2zu(k) *          &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           )                        *          &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )
                       ENDIF
 
-                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
-                                            ( k1* shf(j,i) + k2 * qsws(j,i) )
                    ENDDO
 
+                ENDDO
+
+                IF ( use_surface_fluxes )  THEN
+
+                   DO  j = nys, nyn
+!
+!--                   Treat horizontal default surfaces, up- and downward-facing
+                      DO  l = 0, 1
+                         surf_s = surf_def_h(l)%start_index(j,i)
+                         surf_e = surf_def_h(l)%end_index(j,i)
+                         DO  m = surf_s, surf_e
+                            k = surf_def_h(l)%k(m)
+
+                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                               k2 = 0.61_wp * pt(k,j,i)
+                            ELSE IF ( cloud_physics )  THEN
+                               IF ( ql(k,j,i) == 0.0_wp )  THEN
+                                  k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                                  k2 = 0.61_wp * pt(k,j,i)
+                               ELSE
+                                  theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                                  temp  = theta * t_d_pt(k)
+                                  k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *         &
+                                             ( q(k,j,i) - ql(k,j,i) ) *        &
+                                    ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /   &
+                                    ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *     &
+                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                               ENDIF
+                            ELSE IF ( cloud_droplets )  THEN
+                               k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                               k2 = 0.61_wp * pt(k,j,i)
+                            ENDIF
+
+                            tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *       &
+                                               ( k1 * surf_def_h(l)%shf(m) +   &
+                                                 k2 * surf_def_h(l)%qsws(m) )
+                         ENDDO
+                      ENDDO
+!
+!--                   Treat horizontal natural surfaces
+                      surf_s = surf_lsm_h%start_index(j,i)
+                      surf_e = surf_lsm_h%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_lsm_h%k(m)
+
+                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
+                         ELSE IF ( cloud_physics )  THEN
+                            IF ( ql(k,j,i) == 0.0_wp )  THEN
+                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                               k2 = 0.61_wp * pt(k,j,i)
+                            ELSE
+                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                               temp  = theta * t_d_pt(k)
+                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
+                                             ( q(k,j,i) - ql(k,j,i) ) *        &
+                                    ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /   &
+                                    ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *     &
+                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                            ENDIF
+                         ELSE IF ( cloud_droplets )  THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
+                         ENDIF
+
+                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
+                                               ( k1 * surf_lsm_h%shf(m) +      &
+                                                 k2 * surf_lsm_h%qsws(m) )
+                      ENDDO
+!
+!--                   Treat horizontal urban surfaces
+                      surf_s = surf_usm_h%start_index(j,i)
+                      surf_e = surf_usm_h%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_lsm_h%k(m)
+
+                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
+                         ELSE IF ( cloud_physics )  THEN
+                            IF ( ql(k,j,i) == 0.0_wp )  THEN
+                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                               k2 = 0.61_wp * pt(k,j,i)
+                            ELSE
+                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                               temp  = theta * t_d_pt(k)
+                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
+                                             ( q(k,j,i) - ql(k,j,i) ) *        &
+                                    ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /   &
+                                    ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *     &
+                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                            ENDIF
+                         ELSE IF ( cloud_droplets )  THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
+                         ENDIF
+
+                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
+                                               ( k1 * surf_usm_h%shf(m) +      &
+                                                 k2 * surf_usm_h%qsws(m) )
+                      ENDDO
+
+                   ENDDO
+
                 ENDIF
 
@@ -700,5 +829,608 @@
                    DO  j = nys, nyn
 
-                      k = nzt
+                      surf_s = surf_def_h(2)%start_index(j,i)
+                      surf_e = surf_def_h(2)%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_def_h(2)%k(m)
+
+                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
+                         ELSE IF ( cloud_physics )  THEN
+                            IF ( ql(k,j,i) == 0.0_wp )  THEN
+                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                               k2 = 0.61_wp * pt(k,j,i)
+                            ELSE
+                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                               temp  = theta * t_d_pt(k)
+                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
+                                          ( q(k,j,i) - ql(k,j,i) ) *           &
+                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
+                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
+                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                            ENDIF
+                         ELSE IF ( cloud_droplets )  THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
+                         ENDIF
+
+                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
+                                               ( k1 * surf_def_h(2)%shf(m) +   &
+                                                 k2 * surf_def_h(2)%qsws(m) )
+
+                      ENDDO
+
+                   ENDDO
+
+                ENDIF
+
+             ENDIF
+
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE production_e
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE production_e_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w
+
+       USE cloud_parameters,                                                   &
+           ONLY:  l_d_cp, l_d_r, pt_d_t, t_d_pt
+
+       USE control_parameters,                                                 &
+           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
+                  humidity, kappa, neutral, ocean, pt_reference,               &
+                  rho_reference, use_single_reference_value,                   &
+                  use_surface_fluxes, use_top_fluxes
+
+       USE grid_variables,                                                     &
+           ONLY:  ddx, dx, ddy, dy
+
+       USE indices,                                                            &
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzb, nzt, wall_flags_0
+
+       USE surface_mod,                                                        &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
+                   surf_usm_v 
+
+       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 for different surface type orientation
+       INTEGER(iwp) ::  m       !< running index surface elements
+       INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
+       INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
+
+       REAL(wp)     ::  def         !<
+       REAL(wp)     ::  flag        !< flag to mask topography
+       REAL(wp)     ::  k1          !<
+       REAL(wp)     ::  k2          !<
+       REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
+       REAL(wp)     ::  theta       !<
+       REAL(wp)     ::  temp        !<
+       REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation 
+       REAL(wp)     ::  usvs        !< momentum flux u"v"
+       REAL(wp)     ::  vsus        !< momentum flux v"u"
+       REAL(wp)     ::  wsus        !< momentum flux w"u"
+       REAL(wp)     ::  wsvs        !< momentum flux w"v"
+
+
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudx        !< Gradient of u-component in x-direction
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudy        !< Gradient of u-component in y-direction
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudz        !< Gradient of u-component in z-direction
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdx        !< Gradient of v-component in x-direction
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdy        !< Gradient of v-component in y-direction
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdz        !< Gradient of v-component in z-direction
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdx        !< Gradient of w-component in x-direction
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdy        !< Gradient of w-component in y-direction
+       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdz        !< Gradient of w-component in z-direction
+
+
+       IF ( constant_flux_layer )  THEN
+!
+!--       Calculate TKE production by shear. Calculate gradients at all grid 
+!--       points first, gradients at surface-bounded grid points will be 
+!--       overwritten further below. 
+          DO  k = nzb+1, nzt
+
+             dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+             dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
+                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 
+             dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
+                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 
+
+             dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
+                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx      
+             dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
+                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 
+
+             dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
+                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx      
+             dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
+                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy      
+             dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+          ENDDO
+!
+!--       Compute gradients at north- and south-facing surfaces.
+!--       Note, no vertical natural surfaces so far.
+          DO  l = 0, 1
+!
+!--          Default surfaces
+             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)
+                usvs        = surf_def_v(l)%mom_flux_tke(0,m)
+                wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
+
+                km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp            &
+                                * 0.5_wp * dy
+!
+!--             -1.0 for right-facing wall, 1.0 for left-facing wall
+                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
+                                  BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+                dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
+                dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
+             ENDDO   
+!
+!--          Natural surfaces
+             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)
+                usvs        = surf_lsm_v(l)%mom_flux_tke(0,m)
+                wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
+
+                km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp            &
+                                * 0.5_wp * dy
+!
+!--             -1.0 for right-facing wall, 1.0 for left-facing wall
+                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
+                                  BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+                dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
+                dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
+             ENDDO 
+!
+!--          Urban surfaces
+             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)
+                usvs        = surf_usm_v(l)%mom_flux_tke(0,m)
+                wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
+
+                km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp            &
+                                * 0.5_wp * dy
+!
+!--             -1.0 for right-facing wall, 1.0 for left-facing wall
+                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
+                                  BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+                dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
+                dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
+             ENDDO 
+          ENDDO
+!
+!--       Compute gradients at east- and west-facing walls
+          DO  l = 2, 3
+!
+!--          Default surfaces
+             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)
+                vsus        = surf_def_v(l)%mom_flux_tke(0,m)
+                wsus        = surf_def_v(l)%mom_flux_tke(1,m)
+
+                km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp            &
+                                   * 0.5_wp * dx
+!
+!--             -1.0 for right-facing wall, 1.0 for left-facing wall
+                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
+                                  BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+                dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
+                dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
+             ENDDO 
+!
+!--          Natural surfaces
+             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)
+                vsus        = surf_lsm_v(l)%mom_flux_tke(0,m)
+                wsus        = surf_lsm_v(l)%mom_flux_tke(1,m)
+
+                km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp            &
+                                   * 0.5_wp * dx
+!
+!--             -1.0 for right-facing wall, 1.0 for left-facing wall
+                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
+                                  BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+                dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
+                dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
+             ENDDO 
+!
+!--          Urban surfaces
+             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)
+                vsus        = surf_usm_v(l)%mom_flux_tke(0,m)
+                wsus        = surf_usm_v(l)%mom_flux_tke(1,m)
+
+                km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp            &
+                                   * 0.5_wp * dx
+!
+!--             -1.0 for right-facing wall, 1.0 for left-facing wall
+                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
+                                  BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+                dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
+                dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
+             ENDDO 
+          ENDDO
+!
+!--       Compute gradients at upward-facing walls, first for 
+!--       non-natural default surfaces
+          surf_s = surf_def_h(0)%start_index(j,i)
+          surf_e = surf_def_h(0)%end_index(j,i)
+          DO  m = surf_s, surf_e
+             k = surf_def_h(0)%k(m)
+!
+!--          Please note, actually, an interpolation of u_0 and v_0 
+!--          onto the grid center would be required. However, this 
+!--          would require several data transfers between 2D-grid and
+!--          wall type. The effect of this missing interpolation is 
+!--          negligible. (See also production_e_init).
+             dudz(k)     = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)   
+             dvdz(k)     = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
+
+          ENDDO
+!
+!--       Natural surfaces
+          surf_s = surf_lsm_h%start_index(j,i)
+          surf_e = surf_lsm_h%end_index(j,i)
+          DO  m = surf_s, surf_e
+             k = surf_lsm_h%k(m)
+!
+!--          Please note, actually, an interpolation of u_0 and v_0 
+!--          onto the grid center would be required. However, this 
+!--          would require several data transfers between 2D-grid and
+!--          wall type. The effect of this missing interpolation is 
+!--          negligible. (See also production_e_init).
+             dudz(k)     = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
+             dvdz(k)     = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
+          ENDDO
+!
+!--       Urban surfaces
+          surf_s = surf_usm_h%start_index(j,i)
+          surf_e = surf_usm_h%end_index(j,i)
+          DO  m = surf_s, surf_e
+             k = surf_usm_h%k(m)
+!
+!--          Please note, actually, an interpolation of u_0 and v_0 
+!--          onto the grid center would be required. However, this 
+!--          would require several data transfers between 2D-grid and
+!--          wall type. The effect of this missing interpolation is 
+!--          negligible. (See also production_e_init).
+             dudz(k)     = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
+             dvdz(k)     = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
+          ENDDO
+!
+!--       Compute gradients at downward-facing walls, only for 
+!--       non-natural default surfaces
+          surf_s = surf_def_h(1)%start_index(j,i)
+          surf_e = surf_def_h(1)%end_index(j,i)
+          DO  m = surf_s, surf_e
+             k = surf_def_h(1)%k(m)
+!
+!--          Please note, actually, an interpolation of u_0 and v_0 
+!--          onto the grid center would be required. However, this 
+!--          would require several data transfers between 2D-grid and
+!--          wall type. The effect of this missing interpolation is 
+!--          negligible. (See also production_e_init).
+             dudz(k)     = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
+             dvdz(k)     = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 
+
+          ENDDO
+
+          DO  k = nzb+1, nzt
+
+             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
+                              dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
+                              dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
+                   2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
+
+             IF ( def < 0.0_wp )  def = 0.0_wp
+
+             flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
+
+             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
+
+          ENDDO
+
+       ELSEIF ( use_surface_fluxes )  THEN
+!
+!--       Calculate TKE production by shear. Here, no additional
+!--       wall-bounded code is considered. 
+!--       Why?
+          DO  k = nzb+1, nzt
+
+             dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+             dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
+                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+             dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
+                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+             dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
+                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+             dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
+                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+             dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
+                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+             dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
+                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+             dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
+                              dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
+                              dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
+                   2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
+
+             IF ( def < 0.0_wp )  def = 0.0_wp
+
+             flag        = MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,i), 29 ) ) 
+             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
+
+          ENDDO
+
+       ENDIF
+
+!
+!--    If required, calculate TKE production by buoyancy
+       IF ( .NOT. neutral )  THEN
+
+          IF ( .NOT. humidity )  THEN
+
+             IF ( use_single_reference_value )  THEN
+
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux in
+!--                the bottom and top surface layer
+                   DO  k = nzb+1, nzt
+                      tend(k,j,i) = tend(k,j,i) +                              &
+                                    kh(k,j,i) * g / rho_reference *            &
+                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) *        &
+                                    dd2zu(k) *                                 &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                           )                                   
+                   ENDDO
+
+                ELSE
+
+                   DO  k = nzb+1, nzt
+!
+!--                   Flag 9 is used to mask top fluxes, flag 30 to mask
+!--                   surface fluxes
+                      tend(k,j,i) = tend(k,j,i) -                              &
+                                    kh(k,j,i) * g / pt_reference *             &
+                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           ) *                                 &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )
+
+                   ENDDO
+
+                   IF ( use_surface_fluxes )  THEN
+!
+!--                   Default surfaces, up- and downward-facing
+                      DO  l = 0, 1
+                         surf_s = surf_def_h(l)%start_index(j,i)
+                         surf_e = surf_def_h(l)%end_index(j,i)
+                         DO  m = surf_s, surf_e
+                            k = surf_def_h(l)%k(m)
+                            tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
+                                                           surf_def_h(l)%shf(m) 
+                         ENDDO
+                      ENDDO
+!
+!--                   Natural surfaces
+                      surf_s = surf_lsm_h%start_index(j,i)
+                      surf_e = surf_lsm_h%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_lsm_h%k(m)
+                         tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
+                                                           surf_lsm_h%shf(m) 
+                      ENDDO
+!
+!--                   Urban surfaces
+                      surf_s = surf_usm_h%start_index(j,i)
+                      surf_e = surf_usm_h%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_usm_h%k(m)
+                         tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
+                                                           surf_usm_h%shf(m) 
+                      ENDDO
+                   ENDIF
+
+                   IF ( use_top_fluxes )  THEN
+                      surf_s = surf_def_h(2)%start_index(j,i)
+                      surf_e = surf_def_h(2)%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_def_h(2)%k(m)
+                         tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
+                                                           surf_def_h(2)%shf(m) 
+                      ENDDO
+                   ENDIF
+
+                ENDIF
+
+             ELSE
+
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux in
+!--                the bottom and top surface layer
+                   DO  k = nzb+1, nzt
+                      tend(k,j,i) = tend(k,j,i) +                              &
+                                    kh(k,j,i) * g / prho(k,j,i)       *        &
+                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) *        &
+                                    dd2zu(k) *                                 &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                           )
+                   ENDDO
+
+                ELSE
+
+                   DO  k = nzb+1, nzt
+!
+!--                   Flag 9 is used to mask top fluxes, flag 30 to mask
+!--                   surface fluxes
+                      tend(k,j,i) = tend(k,j,i) -                              &
+                                    kh(k,j,i) * g / pt(k,j,i) *                &
+                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           ) *                                 &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )
+                   ENDDO
+
+                   IF ( use_surface_fluxes )  THEN
+!
+!--                   Default surfaces, up- and downward-facing
+                      DO  l = 0, 1
+                         surf_s = surf_def_h(l)%start_index(j,i)
+                         surf_e = surf_def_h(l)%end_index(j,i)
+                         DO  m = surf_s, surf_e
+                            k = surf_def_h(l)%k(m)
+                            tend(k,j,i) = tend(k,j,i) + g / pt_reference       &
+                                                   * surf_def_h(l)%shf(m)   
+                         ENDDO  
+                      ENDDO
+!
+!--                   Natural surfaces
+                      surf_s = surf_lsm_h%start_index(j,i)
+                      surf_e = surf_lsm_h%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_lsm_h%k(m)
+                         tend(k,j,i) = tend(k,j,i) + g / pt_reference          &
+                                                     * surf_lsm_h%shf(m)   
+                      ENDDO  
+!
+!--                   Urban surfaces
+                      surf_s = surf_usm_h%start_index(j,i)
+                      surf_e = surf_usm_h%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_usm_h%k(m)
+                         tend(k,j,i) = tend(k,j,i) + g / pt_reference          &
+                                                     * surf_usm_h%shf(m)   
+                      ENDDO
+                   ENDIF
+
+                   IF ( use_top_fluxes )  THEN
+                      surf_s = surf_def_h(2)%start_index(j,i)
+                      surf_e = surf_def_h(2)%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_def_h(2)%k(m)
+                         tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
+                                                           surf_def_h(2)%shf(m) 
+                      ENDDO
+                   ENDIF
+
+                ENDIF
+
+             ENDIF
+
+          ELSE
+
+             DO  k = nzb+1, nzt
+!
+!--             Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes
+                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
+                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
+                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
+                                         ) * dd2zu(k) *                        &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           )          *                        &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )
+
+                ELSE IF ( cloud_physics )  THEN
+                   IF ( ql(k,j,i) == 0.0_wp )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
+                   ELSE
+                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                      temp  = theta * t_d_pt(k)
+                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                     &
+                                    ( q(k,j,i) - ql(k,j,i) ) *                 &
+                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /            &
+                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *              &
+                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                   ENDIF
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
+                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
+                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
+                                         ) * dd2zu(k) *                        &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           )          *                        &
+                                      MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )
+                ELSE IF ( cloud_droplets )  THEN
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
+                                     ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +      &
+                                       k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -      &
+                                       pt(k,j,i) * ( ql(k+1,j,i) -             &
+                                                     ql(k-1,j,i) ) ) * dd2zu(k)&
+                                    * MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 30 )  &
+                                           )                                   &
+                                    * MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 9 )   &
+                                           )    
+                ENDIF
+             ENDDO
+
+             IF ( use_surface_fluxes )  THEN
+!
+!--             Treat horizontal default surfaces, up- and downward-facing
+                DO  l = 0, 1
+                   surf_s = surf_def_h(l)%start_index(j,i)
+                   surf_e = surf_def_h(l)%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k = surf_def_h(l)%k(m)
 
                       IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
@@ -710,12 +1442,12 @@
                             k2 = 0.61_wp * pt(k,j,i)
                          ELSE
-                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                            temp  = theta * t_d_pt(k)
-                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
-                                          ( q(k,j,i) - ql(k,j,i) ) *           &
-                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
-                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
-                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                           theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                           temp  = theta * t_d_pt(k)
+                           k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                &
+                                      ( q(k,j,i) - ql(k,j,i) ) *               &
+                             ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /          &
+                             ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *            & 
+                             ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                           k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                       ELSE IF ( cloud_droplets )  THEN
@@ -724,536 +1456,113 @@
                       ENDIF
 
-                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
-                                            ( k1* tswst(j,i) + k2 * qswst(j,i) )
+                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *             &
+                                         ( k1 * surf_def_h(l)%shf(m) +         &
+                                           k2 * surf_def_h(l)%qsws(m) )
                    ENDDO
-
-                ENDIF
+                ENDDO
+!
+!--             Treat horizontal natural surfaces
+                surf_s = surf_lsm_h%start_index(j,i)
+                surf_e = surf_lsm_h%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_lsm_h%k(m)
+
+                   IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                       k2 = 0.61_wp * pt(k,j,i)
+                   ELSE IF ( cloud_physics )  THEN
+                       IF ( ql(k,j,i) == 0.0_wp )  THEN
+                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                          k2 = 0.61_wp * pt(k,j,i)
+                       ELSE
+                          theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                          temp  = theta * t_d_pt(k)
+                          k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
+                                        ( q(k,j,i) - ql(k,j,i) ) *             &
+                               ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+                               ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
+                               ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                          k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                      ENDIF
+                   ELSE IF ( cloud_droplets )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
+                   ENDIF
+
+                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
+                                            ( k1 * surf_lsm_h%shf(m) +         &
+                                              k2 * surf_lsm_h%qsws(m) )
+                ENDDO
+!
+!--             Treat horizontal urban surfaces
+                surf_s = surf_usm_h%start_index(j,i)
+                surf_e = surf_usm_h%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_usm_h%k(m)
+
+                   IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                       k2 = 0.61_wp * pt(k,j,i)
+                   ELSE IF ( cloud_physics )  THEN
+                       IF ( ql(k,j,i) == 0.0_wp )  THEN
+                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                          k2 = 0.61_wp * pt(k,j,i)
+                       ELSE
+                          theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                          temp  = theta * t_d_pt(k)
+                          k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
+                                        ( q(k,j,i) - ql(k,j,i) ) *             &
+                               ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+                               ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
+                               ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                          k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                      ENDIF
+                   ELSE IF ( cloud_droplets )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
+                   ENDIF
+
+                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
+                                            ( k1 * surf_usm_h%shf(m) +         &
+                                              k2 * surf_usm_h%qsws(m) )
+                ENDDO
 
              ENDIF
 
-          ENDIF
-
-       ENDDO
-
-    END SUBROUTINE production_e
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE production_e_ij( i, j )
-
-       USE arrays_3d,                                                          &
-           ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, qsws, qswst, shf,      &
-                  tend, tswst, u, v, vpt, w
-
-       USE cloud_parameters,                                                   &
-           ONLY:  l_d_cp, l_d_r, pt_d_t, t_d_pt
-
-       USE control_parameters,                                                 &
-           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
-                  humidity, kappa, neutral, ocean, pt_reference,               &
-                  rho_reference, use_single_reference_value,                   &
-                  use_surface_fluxes, use_top_fluxes
-
-       USE grid_variables,                                                     &
-           ONLY:  ddx, dx, ddy, dy, wall_e_x, wall_e_y
-
-       USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner,                   &
-                  nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i           !<
-       INTEGER(iwp) ::  j           !<
-       INTEGER(iwp) ::  k           !<
-
-       REAL(wp)     ::  def         !<
-       REAL(wp)     ::  dudx        !<
-       REAL(wp)     ::  dudy        !<
-       REAL(wp)     ::  dudz        !<
-       REAL(wp)     ::  dvdx        !<
-       REAL(wp)     ::  dvdy        !<
-       REAL(wp)     ::  dvdz        !<
-       REAL(wp)     ::  dwdx        !<
-       REAL(wp)     ::  dwdy        !<
-       REAL(wp)     ::  dwdz        !<
-       REAL(wp)     ::  k1          !<
-       REAL(wp)     ::  k2          !<
-       REAL(wp)     ::  km_neutral  !<
-       REAL(wp)     ::  theta       !<
-       REAL(wp)     ::  temp        !<
-
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  usvs  !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus  !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs  !<
-
-!
-!--    Calculate TKE production by shear
-       DO  k = nzb_diff_s_outer(j,i), nzt
-
-          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
-                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2    &
-                + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-          IF ( def < 0.0_wp )  def = 0.0_wp
-
-          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-       ENDDO
-
-       IF ( constant_flux_layer )  THEN
-
-          IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) )  THEN
-
-!
-!--          Position beneath wall
-!--          (2) - Will allways be executed.
-!--          'bottom and wall: use u_0,v_0 and wall functions'
-             k = nzb_diff_s_inner(j,i)-1
-
-             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-             dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                               u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-             dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                               v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
-
-             IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
-!
-!--             Inconsistency removed: as the thermal stratification 
-!--             is not taken into account for the evaluation of the 
-!--             wall fluxes at vertical walls, the eddy viscosity km 
-!--             must not be used for the evaluation of the velocity 
-!--             gradients dudy and dwdy
-!--             Note: The validity of the new method has not yet 
-!--                   been shown, as so far no suitable data for a 
-!--                   validation has been available
-                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
-                                    usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
-                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
-                                    wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
-                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
-                             0.5_wp * dy
-                IF ( km_neutral > 0.0_wp )  THEN
-                   dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
-                   dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
-                ELSE
-                   dudy = 0.0_wp
-                   dwdy = 0.0_wp
-                ENDIF
-             ELSE
-                dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-             ENDIF
-
-             IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
-!
-!--             Inconsistency removed: as the thermal stratification 
-!--             is not taken into account for the evaluation of the 
-!--             wall fluxes at vertical walls, the eddy viscosity km 
-!--             must not be used for the evaluation of the velocity 
-!--             gradients dvdx and dwdx
-!--             Note: The validity of the new method has not yet 
-!--                   been shown, as so far no suitable data for a 
-!--                   validation has been available
-                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
-                                    vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
-                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
-                                    wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
-                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 
-                             0.5_wp * dx
-                IF ( km_neutral > 0.0_wp )  THEN
-                   dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
-                   dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
-                ELSE
-                   dvdx = 0.0_wp
-                   dwdx = 0.0_wp
-                ENDIF
-             ELSE
-                dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-             ENDIF
-
-             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                   dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-             IF ( def < 0.0_wp )  def = 0.0_wp
-
-             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-!
-!--          (3) - will be executed only, if there is at least one level
-!--          between (2) and (4), i.e. the topography must have a
-!--          minimum height of 2 dz. Wall fluxes for this case have
-!--          already been calculated for (2).
-!--          'wall only: use wall functions'
-             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
-
-                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-                dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
-
-                IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
-!
-!--                Inconsistency removed: as the thermal stratification 
-!--                is not taken into account for the evaluation of the 
-!--                wall fluxes at vertical walls, the eddy viscosity km 
-!--                must not be used for the evaluation of the velocity 
-!--                gradients dudy and dwdy
-!--                Note: The validity of the new method has not yet 
-!--                      been shown, as so far no suitable data for a 
-!--                      validation has been available
-                   km_neutral = kappa * ( usvs(k)**2 + & 
-                                          wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
-                   IF ( km_neutral > 0.0_wp )  THEN
-                      dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
-                      dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
-                   ELSE
-                      dudy = 0.0_wp
-                      dwdy = 0.0_wp
-                   ENDIF
-                ELSE
-                   dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                ENDIF
-
-                IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
-!
-!--                Inconsistency removed: as the thermal stratification 
-!--                is not taken into account for the evaluation of the 
-!--                wall fluxes at vertical walls, the eddy viscosity km 
-!--                must not be used for the evaluation of the velocity 
-!--                gradients dvdx and dwdx
-!--                Note: The validity of the new method has not yet 
-!--                      been shown, as so far no suitable data for a 
-!--                      validation has been available
-                   km_neutral = kappa * ( vsus(k)**2 + & 
-                                          wsus(k)**2 )**0.25_wp * 0.5_wp * dx
-                   IF ( km_neutral > 0.0_wp )  THEN
-                      dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
-                      dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
-                   ELSE
-                      dvdx = 0.0_wp
-                      dwdx = 0.0_wp
-                   ENDIF
-                ELSE
-                   dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                ENDIF
-
-                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0_wp )  def = 0.0_wp
-
-                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-             ENDDO
-
-!
-!--          (4) - will allways be executed.
-!--          'special case: free atmosphere' (as for case (0))
-             k = nzb_diff_s_outer(j,i)-1
-
-             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +        &
-                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                   dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-             IF ( def < 0.0_wp )  def = 0.0_wp
-
-             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-          ELSE
-
-!
-!--          Position without adjacent wall
-!--          (1) - will allways be executed. 
-!--          'bottom only: use u_0,v_0'
-             k = nzb_diff_s_inner(j,i)-1
-
-             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                 u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-
-             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                 v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
-                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
-                   + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-             IF ( def < 0.0_wp )  def = 0.0_wp
-
-             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-          ENDIF
-
-       ELSEIF ( use_surface_fluxes )  THEN
-
-          k = nzb_diff_s_outer(j,i)-1
-
-          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-          IF ( def < 0.0_wp )  def = 0.0_wp
-
-          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-       ENDIF
-
-!
-!--    If required, calculate TKE production by buoyancy
-       IF ( .NOT. neutral )  THEN
-
-          IF ( .NOT. humidity )  THEN
-
-             IF ( use_single_reference_value )  THEN
-
-                IF ( ocean )  THEN
-!
-!--                So far in the ocean no special treatment of density flux in
-!--                the bottom and top surface layer
-                   DO  k = nzb_s_inner(j,i)+1, nzt
-                      tend(k,j,i) = tend(k,j,i) +                   &
-                                    kh(k,j,i) * g / rho_reference * &
-                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
-                   ENDDO
-
-                ELSE
-
-                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                      tend(k,j,i) = tend(k,j,i) -                  &
-                                    kh(k,j,i) * g / pt_reference * &
-                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
-                   ENDDO
-
-                   IF ( use_surface_fluxes )  THEN
-                      k = nzb_diff_s_inner(j,i)-1
-                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
-                   ENDIF
-
-                   IF ( use_top_fluxes )  THEN
-                      k = nzt
-                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
-                   ENDIF
-
-                ENDIF
-
-             ELSE
-
-                IF ( ocean )  THEN
-!
-!--                So far in the ocean no special treatment of density flux in
-!--                the bottom and top surface layer
-                   DO  k = nzb_s_inner(j,i)+1, nzt
-                      tend(k,j,i) = tend(k,j,i) +                              &
-                                    kh(k,j,i) * g / prho(k,j,i) *              &
-                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
-                   ENDDO
-
-                ELSE
-
-                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                      tend(k,j,i) = tend(k,j,i) -               &
-                                    kh(k,j,i) * g / pt(k,j,i) * &
-                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
-                   ENDDO
-
-                   IF ( use_surface_fluxes )  THEN
-                      k = nzb_diff_s_inner(j,i)-1
-                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
-                   ENDIF
-
-                   IF ( use_top_fluxes )  THEN
-                      k = nzt
-                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
-                   ENDIF
-
-                ENDIF
-
-             ENDIF
-
-          ELSE
-
-             DO  k = nzb_diff_s_inner(j,i), nzt_diff
-
-                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                   k2 = 0.61_wp * pt(k,j,i)
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
-                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
-                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
-                                         ) * dd2zu(k)
-                ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0_wp )  THEN
+             IF ( use_top_fluxes )  THEN
+                surf_s = surf_def_h(2)%start_index(j,i)
+                surf_e = surf_def_h(2)%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_def_h(2)%k(m)
+
+
+
+                   IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                       k2 = 0.61_wp * pt(k,j,i)
-                   ELSE
-                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                      temp  = theta * t_d_pt(k)
-                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
-                                    ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
-                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
+                   ELSE IF ( cloud_physics )  THEN
+                      IF ( ql(k,j,i) == 0.0_wp )  THEN
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
+                      ELSE
+                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                         temp  = theta * t_d_pt(k)
+                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
+                                    ( q(k,j,i) - ql(k,j,i) ) *                 &
+                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /            &
+                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *              &
                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                         k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                      ENDIF
+                   ELSE IF ( cloud_droplets )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
                    ENDIF
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
-                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
-                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
-                                         ) * dd2zu(k)
-                ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61_wp * pt(k,j,i)
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *  &
-                                     ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +    &
-                                       k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -    &
-                                       pt(k,j,i) * ( ql(k+1,j,i) -           &
-                                                     ql(k-1,j,i) ) ) * dd2zu(k)
-                ENDIF
-             ENDDO
-
-             IF ( use_surface_fluxes )  THEN
-                k = nzb_diff_s_inner(j,i)-1
-
-                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                   k2 = 0.61_wp * pt(k,j,i)
-                ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0_wp )  THEN
-                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                      k2 = 0.61_wp * pt(k,j,i)
-                   ELSE
-                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                      temp  = theta * t_d_pt(k)
-                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
-                                    ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
-                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
-                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
-                   ENDIF
-                ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61_wp * pt(k,j,i)
-                ENDIF
-
-                tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
-                                            ( k1* shf(j,i) + k2 * qsws(j,i) )
-             ENDIF
-
-             IF ( use_top_fluxes )  THEN
-                k = nzt
-
-                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                   k2 = 0.61_wp * pt(k,j,i)
-                ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0_wp )  THEN
-                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                      k2 = 0.61_wp * pt(k,j,i)
-                   ELSE
-                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                      temp  = theta * t_d_pt(k)
-                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
-                                    ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
-                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
-                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
-                   ENDIF
-                ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61_wp * pt(k,j,i)
-                ENDIF
-
-                tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
-                                            ( k1* tswst(j,i) + k2 * qswst(j,i) )
+
+                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
+                               ( k1* surf_def_h(2)%shf(m) +                    &
+                                 k2 * surf_def_h(2)%qsws(m) )
+                ENDDO
+
              ENDIF
 
@@ -1273,5 +1582,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  kh, km, u, us, usws, v, vsws, zu
+           ONLY:  kh, km, u, v, zu
 
        USE control_parameters,                                                 &
@@ -1282,20 +1591,16 @@
                   nzb_v_inner
 
+       USE surface_mod,                                                        &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h 
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i   !<
-       INTEGER(iwp) ::  j   !<
-       INTEGER(iwp) ::  ku  !<
-       INTEGER(iwp) ::  kv  !<
+       INTEGER(iwp) ::  i   !< grid index x-direction
+       INTEGER(iwp) ::  j   !< grid index y-direction
+       INTEGER(iwp) ::  k   !< grid index z-direction
+       INTEGER(iwp) ::  l   !< running index surface type (up- or downward-facing)
+       INTEGER(iwp) ::  m   !< running index surface elements
 
        IF ( constant_flux_layer )  THEN
-
-          IF ( first_call )  THEN
-             ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
-             u_0 = 0.0_wp   ! just to avoid access of uninitialized memory
-             v_0 = 0.0_wp   ! within exchange_horiz_2d
-             first_call = .FALSE.
-          ENDIF
-
 !
 !--       Calculate a virtual velocity at the surface in a way that the
@@ -1306,34 +1611,124 @@
 !--       (otherwise the timestep may be significantly reduced by large
 !--       surface winds).
-!--       Upper bounds are nxr+1 and nyn+1 because otherwise these values are
 !--       not available in case of non-cyclic boundary conditions.
 !--       WARNING: the exact analytical solution would require the determination
 !--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
-          !$OMP PARALLEL DO PRIVATE( ku, kv )
-          DO  i = nxl, nxr+1
-             DO  j = nys, nyn+1
-
-                ku = nzb_u_inner(j,i)+1
-                kv = nzb_v_inner(j,i)+1
-
-                u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
-                                 ( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) +    &
-                                   1.0E-20_wp )
-!                                  ( us(j,i) * kappa * zu(1) )
-                v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
-                                 ( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) +    &
-                                   1.0E-20_wp )
-!                                  ( us(j,i) * kappa * zu(1) )
-
-                IF ( ABS( u(ku+1,j,i) - u_0(j,i) )  > &
-                     ABS( u(ku+1,j,i) - u(ku-1,j,i) ) )  u_0(j,i) = u(ku-1,j,i)
-                IF ( ABS( v(kv+1,j,i) - v_0(j,i) )  > &
-                     ABS( v(kv+1,j,i) - v(kv-1,j,i) ) )  v_0(j,i) = v(kv-1,j,i)
-
-             ENDDO
+!--       Default surfaces, upward-facing
+          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+          DO  m = 1, surf_def_h(0)%ns
+
+             i = surf_def_h(0)%i(m)            
+             j = surf_def_h(0)%j(m)
+             k = surf_def_h(0)%k(m)
+!
+!--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 
+!--          and km are not on the same grid. Actually, a further 
+!--          interpolation of km onto the u/v-grid is necessary. However, the
+!--          effect of this error is negligible. 
+             surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *       &
+                                        ( zu(k+1)    - zu(k-1)    )  /         &
+                                        ( km(k,j,i)  + 1.0E-20_wp )  
+             surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *       &
+                                        ( zu(k+1)    - zu(k-1)    )  /         &
+                                        ( km(k,j,i)  + 1.0E-20_wp )  
+
+             IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                  &
+                  ABS( u(k+1,j,i) - u(k-1,j,i)           )                     &
+                )  surf_def_h(0)%u_0(m) = u(k-1,j,i)
+
+             IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) )  >                  &
+                  ABS( v(k+1,j,i) - v(k-1,j,i)           )                     &
+                )  surf_def_h(0)%v_0(m) = v(k-1,j,i)
+
           ENDDO
-
-          CALL exchange_horiz_2d( u_0 )
-          CALL exchange_horiz_2d( v_0 )
+!
+!--       Default surfaces, downward-facing
+          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+          DO  m = 1, surf_def_h(1)%ns
+
+             i = surf_def_h(1)%i(m)            
+             j = surf_def_h(1)%j(m)
+             k = surf_def_h(1)%k(m)
+!
+!--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 
+!--          and km are not on the same grid. Actually, a further 
+!--          interpolation of km onto the u/v-grid is necessary. However, the
+!--          effect of this error is negligible. 
+!--          In case of downward-facing surfaces, gradient is calculated 
+!--          between u_0 and u(k-1).
+             surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *       &
+                                        ( zu(k+1)    - zu(k-1)    )  /         &
+                                        ( km(k,j,i)  + 1.0E-20_wp )  
+             surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *       &
+                                        ( zu(k+1)    - zu(k-1)    )  /         &
+                                        ( km(k,j,i)  + 1.0E-20_wp )  
+
+             IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >                  &
+                  ABS( u(k+1,j,i)           - u(k-1,j,i) )                     &
+                )  surf_def_h(1)%u_0(m) = u(k+1,j,i)
+
+             IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >                  &
+                  ABS( v(k+1,j,i)           - v(k-1,j,i) )                     &
+                )  surf_def_h(1)%v_0(m) = v(k+1,j,i)
+
+          ENDDO
+!
+!--       Natural surfaces, upward-facing
+          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+          DO  m = 1, surf_lsm_h%ns
+
+             i = surf_lsm_h%i(m)            
+             j = surf_lsm_h%j(m)
+             k = surf_lsm_h%k(m)
+!
+!--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 
+!--          and km are not on the same grid. Actually, a further 
+!--          interpolation of km onto the u/v-grid is necessary. However, the
+!--          effect of this error is negligible. 
+             surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)      *        &
+                                           ( zu(k+1)   - zu(k-1)    )  /       &
+                                           ( km(k,j,i) + 1.0E-20_wp )  
+             surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)      *        &
+                                           ( zu(k+1)   - zu(k-1)    )  /       &
+                                           ( km(k,j,i) + 1.0E-20_wp )
+
+             IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                     &
+                  ABS( u(k+1,j,i) - u(k-1,j,i)   )                             &
+                )  surf_lsm_h%u_0(m) = u(k-1,j,i)
+
+             IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) )  >                     &
+                  ABS( v(k+1,j,i) - v(k-1,j,i)   )                             &
+                )  surf_lsm_h%v_0(m) = v(k-1,j,i)
+
+          ENDDO
+!
+!--       Urban surfaces, upward-facing
+          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+          DO  m = 1, surf_usm_h%ns
+
+             i = surf_usm_h%i(m)            
+             j = surf_usm_h%j(m)
+             k = surf_usm_h%k(m)
+!
+!--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 
+!--          and km are not on the same grid. Actually, a further 
+!--          interpolation of km onto the u/v-grid is necessary. However, the
+!--          effect of this error is negligible. 
+             surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)      *        &
+                                           ( zu(k+1)   - zu(k-1)    )  /       &
+                                           ( km(k,j,i) + 1.0E-20_wp )  
+             surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)      *        &
+                                           ( zu(k+1)   - zu(k-1)    )  /       &
+                                           ( km(k,j,i) + 1.0E-20_wp )
+
+             IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                     &
+                  ABS( u(k+1,j,i) - u(k-1,j,i)   )                             &
+                )  surf_usm_h%u_0(m) = u(k-1,j,i)
+
+             IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) )  >                     &
+                  ABS( v(k+1,j,i) - v(k-1,j,i)   )                             &
+                )  surf_usm_h%v_0(m) = v(k-1,j,i)
+
+          ENDDO
 
        ENDIF
