Index: palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_s.f90	(revision 2119)
+++ palm/trunk/SOURCE/diffusion_s.f90	(revision 2232)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Adjustments to new topography and surface concept
 ! 
 ! Former revisions:
@@ -110,5 +110,13 @@
 !> Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux )
+    SUBROUTINE diffusion_s( s, s_flux_def_h_up,    s_flux_def_h_down,          &
+                               s_flux_t,                                       &
+                               s_flux_lsm_h_up,    s_flux_usm_h_up,            &
+                               s_flux_def_v_north, s_flux_def_v_south,         &
+                               s_flux_def_v_east,  s_flux_def_v_west,          &
+                               s_flux_lsm_v_north, s_flux_lsm_v_south,         &
+                               s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
+                               s_flux_usm_v_north, s_flux_usm_v_south,         &
+                               s_flux_usm_v_east,  s_flux_usm_v_west )
 
        USE arrays_3d,                                                          &
@@ -119,19 +127,50 @@
        
        USE grid_variables,                                                     &
-           ONLY:  ddx2, ddy2, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
+           ONLY:  ddx2, ddy2
        
        USE indices,                                                            &
            ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb,             &
-                  nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, nzt, nzt_diff
+                  nzt, wall_flags_0
        
        USE kinds
 
+       USE surface_mod,                                                        &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
+                   surf_usm_v 
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                 !< 
-       INTEGER(iwp) ::  j                 !< 
-       INTEGER(iwp) ::  k                 !< 
-       REAL(wp)     ::  wall_s_flux(0:4)  !< 
-       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t !<
+       INTEGER(iwp) ::  i             !< running index x direction
+       INTEGER(iwp) ::  j             !< running index y direction
+       INTEGER(iwp) ::  k             !< running index z direction
+       INTEGER(iwp) ::  m             !< running index surface elements
+       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
+       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
+
+       REAL(wp) ::  flag              !< flag to mask topography grid points
+       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
+       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point 
+       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
+       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point 
+       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
+       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface  
+
+       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces 
+       REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical natural-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical natural-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
 #if defined( __nopointer )
        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !< 
@@ -144,95 +183,236 @@
 !
 !--          Compute horizontal diffusion
-             DO  k = nzb_s_outer(j,i)+1, nzt
+             DO  k = nzb+1, nzt
+!
+!--             Predetermine flag to mask topography and wall-bounded grid points
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
+!
+!--             Predetermine flag to mask wall-bounded grid points, equivalent to
+!--             former s_outer array
+                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
+                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
+                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
+                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
 
                 tend(k,j,i) = tend(k,j,i)                                      &
                                           + 0.5_wp * (                         &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
-                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                                                     ) * ddx2                  &
+                        mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )               &
+                                   * ( s(k,j,i+1) - s(k,j,i)   )               &
+                      - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )               &
+                                   * ( s(k,j,i)   - s(k,j,i-1) )               &
+                                                     ) * ddx2 * flag           &
                                           + 0.5_wp * (                         &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
-                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                                                     ) * ddy2
-             ENDDO
-
-!
-!--          Apply prescribed horizontal wall heatflux where necessary
-             IF ( ( wall_w_x(j,i) /= 0.0_wp ) .OR. ( wall_w_y(j,i) /= 0.0_wp ) &
-                )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
-
+                        mask_north * ( kh(k,j,i) + kh(k,j+1,i) )               &
+                                   * ( s(k,j+1,i) - s(k,j,i)   )               &
+                      - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )               &
+                                   * ( s(k,j,i)   - s(k,j-1,i) )               &
+                                                     ) * ddy2 * flag
+             ENDDO
+
+!
+!--          Apply prescribed horizontal wall heatflux where necessary. First,
+!--          determine start and end index for respective (j,i)-index. Please
+!--          note, in the flat case following loop will not be entered, as
+!--          surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
+!--          so far. 
+!--          First, for default-type surfaces 
+!--          North-facing vertical default-type surfaces
+             surf_s = surf_def_v(0)%start_index(j,i)
+             surf_e = surf_def_v(0)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_def_v(0)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy2
+             ENDDO
+!
+!--          South-facing vertical default-type surfaces
+             surf_s = surf_def_v(1)%start_index(j,i)
+             surf_e = surf_def_v(1)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_def_v(1)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy2
+             ENDDO
+!
+!--          East-facing vertical default-type surfaces
+             surf_s = surf_def_v(2)%start_index(j,i)
+             surf_e = surf_def_v(2)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_def_v(2)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx2
+             ENDDO
+!
+!--          West-facing vertical default-type surfaces
+             surf_s = surf_def_v(3)%start_index(j,i)
+             surf_e = surf_def_v(3)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_def_v(3)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx2
+             ENDDO
+!
+!--          Now, for natural-type surfaces.
+!--          North-facing
+             surf_s = surf_lsm_v(0)%start_index(j,i)
+             surf_e = surf_lsm_v(0)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_lsm_v(0)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy2
+             ENDDO
+!
+!--          South-facing
+             surf_s = surf_lsm_v(1)%start_index(j,i)
+             surf_e = surf_lsm_v(1)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_lsm_v(1)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy2
+             ENDDO
+!
+!--          East-facing
+             surf_s = surf_lsm_v(2)%start_index(j,i)
+             surf_e = surf_lsm_v(2)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_lsm_v(2)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx2
+             ENDDO
+!
+!--          West-facing
+             surf_s = surf_lsm_v(3)%start_index(j,i)
+             surf_e = surf_lsm_v(3)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_lsm_v(3)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx2
+             ENDDO
+!
+!--          Now, for urban-type surfaces.
+!--          North-facing
+             surf_s = surf_usm_v(0)%start_index(j,i)
+             surf_e = surf_usm_v(0)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_usm_v(0)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy2
+             ENDDO
+!
+!--          South-facing
+             surf_s = surf_usm_v(1)%start_index(j,i)
+             surf_e = surf_usm_v(1)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_usm_v(1)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy2
+             ENDDO
+!
+!--          East-facing
+             surf_s = surf_usm_v(2)%start_index(j,i)
+             surf_e = surf_usm_v(2)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_usm_v(2)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx2
+             ENDDO
+!
+!--          West-facing
+             surf_s = surf_usm_v(3)%start_index(j,i)
+             surf_e = surf_usm_v(3)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k           = surf_usm_v(3)%k(m)
+                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx2
+             ENDDO
+
+!
+!--          Compute vertical diffusion. In case that surface fluxes have been
+!--          prescribed or computed at bottom and/or top, index k starts/ends at
+!--          nzb+2 or nzt-1, respectively. Model top is also mask if top flux
+!--          is given.
+             DO  k = nzb+1, nzt
+!
+!--             Determine flags to mask topography below and above. Flag 0 is 
+!--             used to mask topography in general, and flag 8 implies 
+!--             information about use_surface_fluxes. Flag 9 is used to control 
+!--             flux at model top. 
+                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
+                                     BTEST( wall_flags_0(k-1,j,i), 8 ) ) 
+                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
+                                     BTEST( wall_flags_0(k+1,j,i), 8 ) ) *     &
+                              MERGE( 1.0_wp, 0.0_wp,                           &
+                                     BTEST( wall_flags_0(k+1,j,i), 9 ) ) 
+                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
+                                     BTEST( wall_flags_0(k,j,i), 0 ) )
+
+                tend(k,j,i) = tend(k,j,i)                                      &
+                                       + 0.5_wp * (                            &
+                                      ( kh(k,j,i) + kh(k+1,j,i) ) *            &
+                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
+                                                            * rho_air_zw(k)    &
+                                                            * mask_top         &
+                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
+                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
+                                                            * rho_air_zw(k-1)  &
+                                                            * mask_bottom      &
+                                                  ) * ddzw(k) * drho_air(k)    &
+                                                              * flag
+             ENDDO
+
+!
+!--          Vertical diffusion at horizontal walls.
+             IF ( use_surface_fluxes )  THEN
+!
+!--             Default-type surfaces, upward-facing               
+                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)
+                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)              &
+                                       * ddzw(k) * drho_air(k)
+
+                ENDDO
+!
+!--             Default-type surfaces, downward-facing               
+                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)
+                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)            &
+                                       * ddzw(k) * drho_air(k)
+
+                ENDDO
+!
+!--             Natural-type surfaces, upward-facing  
+                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) + s_flux_lsm_h_up(m)              &
+                                       * ddzw(k) * drho_air(k)
+
+                ENDDO
+!
+!--             Urban-type surfaces, upward-facing     
+                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) + s_flux_usm_h_up(m)              &
+                                       * ddzw(k) * drho_air(k)
+
+                ENDDO
+
+             ENDIF
+!
+!--          Vertical diffusion at the last computational gridpoint along z-direction
+             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)                                   &
-                                                + ( fwxp(j,i) * 0.5_wp *       &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
-                        + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
-                                                   -fwxm(j,i) * 0.5_wp *       &
-                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                        + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
-                                                  ) * ddx2                     &
-                                                + ( fwyp(j,i) * 0.5_wp *       &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
-                        + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
-                                                   -fwym(j,i) * 0.5_wp *       &
-                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                        + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
-                                                  ) * ddy2
+                           + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
                 ENDDO
              ENDIF
 
-!
-!--          Compute vertical diffusion. In case that surface fluxes have been
-!--          prescribed or computed at bottom and/or top, index k starts/ends at
-!--          nzb+2 or nzt-1, respectively.
-             DO  k = nzb_diff_s_inner(j,i), nzt_diff
-
-                tend(k,j,i) = tend(k,j,i)                                      &
-                                       + 0.5_wp * (                            &
-            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
-                                                            * rho_air_zw(k)    &
-          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                                            * rho_air_zw(k-1)  &
-                                                  ) * ddzw(k) * drho_air(k)
-             ENDDO
-
-!
-!--          Vertical diffusion at the first computational gridpoint along
-!--          z-direction
-             IF ( use_surface_fluxes )  THEN
-
-                k = nzb_s_inner(j,i)+1
-
-                tend(k,j,i) = tend(k,j,i)                                      &
-                                       + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )  &
-                                                  * ( s(k+1,j,i)-s(k,j,i) )    &
-                                                  * ddzu(k+1)                  &
-                                                  * rho_air_zw(k)              &
-                                           + s_flux_b(j,i)                     &
-                                         ) * ddzw(k) * drho_air(k)
-
-             ENDIF
-
-!
-!--          Vertical diffusion at the last computational gridpoint along
-!--          z-direction
-             IF ( use_top_fluxes )  THEN
-
-                k = nzt
-
-                tend(k,j,i) = tend(k,j,i)                                      &
-                                       + ( - s_flux_t(j,i)                     &
-                                           - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )&
-                                                    * ( s(k,j,i)-s(k-1,j,i) )  &
-                                                    * ddzu(k)                  &
-                                                    * rho_air_zw(k-1)          &
-                                         ) * ddzw(k) * drho_air(k)
-
-             ENDIF
-
           ENDDO
        ENDDO
 
     END SUBROUTINE diffusion_s
-
 
 !------------------------------------------------------------------------------!
@@ -241,5 +421,14 @@
 !> Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
+    SUBROUTINE diffusion_s_ij( i, j, s,                                        &
+                               s_flux_def_h_up,    s_flux_def_h_down,          &
+                               s_flux_t,                                       &
+                               s_flux_lsm_h_up,    s_flux_usm_h_up,            &
+                               s_flux_def_v_north, s_flux_def_v_south,         &
+                               s_flux_def_v_east,  s_flux_def_v_west,          &
+                               s_flux_lsm_v_north, s_flux_lsm_v_south,         &
+                               s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
+                               s_flux_usm_v_north, s_flux_usm_v_south,         &
+                               s_flux_usm_v_east,  s_flux_usm_v_west )       
 
        USE arrays_3d,                                                          &
@@ -250,20 +439,49 @@
        
        USE grid_variables,                                                     &
-           ONLY:  ddx2, ddy2, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
+           ONLY:  ddx2, ddy2
        
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_diff_s_inner, nzb_s_inner,  &
-                  nzb_s_outer, nzt, nzt_diff
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
        
        USE kinds
 
+       USE surface_mod,                                                        &
+           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
+                   surf_usm_v 
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                 !< 
-       INTEGER(iwp) ::  j                 !< 
-       INTEGER(iwp) ::  k                 !< 
-       REAL(wp)     ::  wall_s_flux(0:4)  !< 
-       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b  !< 
-       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_t  !< 
+       INTEGER(iwp) ::  i             !< running index x direction
+       INTEGER(iwp) ::  j             !< running index y direction
+       INTEGER(iwp) ::  k             !< running index z direction
+       INTEGER(iwp) ::  m             !< running index surface elements
+       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
+       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
+
+       REAL(wp) ::  flag              !< flag to mask topography grid points
+       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
+       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point 
+       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
+       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point 
+       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
+       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface  
+
+       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces 
+       REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
+       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
 #if defined( __nopointer )
        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s !< 
@@ -274,84 +492,232 @@
 !
 !--    Compute horizontal diffusion
-       DO  k = nzb_s_outer(j,i)+1, nzt
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography and wall-bounded grid points
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
+!
+!--       Predetermine flag to mask wall-bounded grid points, equivalent to
+!--       former s_outer array
+          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
+          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
+          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
+          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
+!
+!--       Finally, determine flag to mask both topography itself as well
+!--       as wall-bounded grid points, which will be treated further below
 
           tend(k,j,i) = tend(k,j,i)                                            &
                                           + 0.5_wp * (                         &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
-                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                                                     ) * ddx2                  &
+                            mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )           &
+                                       * ( s(k,j,i+1) - s(k,j,i)   )           &
+                          - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )           &
+                                       * ( s(k,j,i)   - s(k,j,i-1) )           &
+                                                     ) * ddx2 * flag           &
                                           + 0.5_wp * (                         &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
-                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                                                     ) * ddy2
-       ENDDO
-
-!
-!--    Apply prescribed horizontal wall heatflux where necessary
-       IF ( ( wall_w_x(j,i) /= 0.0_wp ) .OR. ( wall_w_y(j,i) /= 0.0_wp ) )     &
-       THEN
-          DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
-
-             tend(k,j,i) = tend(k,j,i)                                         &
-                                                + ( fwxp(j,i) * 0.5_wp *       &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
-                        + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
-                                                   -fwxm(j,i) * 0.5_wp *       &
-                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                        + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
-                                                  ) * ddx2                     &
-                                                + ( fwyp(j,i) * 0.5_wp *       &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
-                        + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
-                                                   -fwym(j,i) * 0.5_wp *       &
-                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                        + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
-                                                  ) * ddy2
+                            mask_north * ( kh(k,j,i) + kh(k,j+1,i) )           &
+                                       * ( s(k,j+1,i) - s(k,j,i)   )           &
+                          - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )           &
+                                       * ( s(k,j,i)  - s(k,j-1,i)  )           &
+                                                     ) * ddy2 * flag
+       ENDDO
+
+!
+!--    Apply prescribed horizontal wall heatflux where necessary. First,
+!--    determine start and end index for respective (j,i)-index. Please
+!--    note, in the flat case following loops will not be entered, as
+!--    surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
+!--    so far. 
+!--    First, for default-type surfaces
+!--    North-facing vertical default-type surfaces
+       surf_s = surf_def_v(0)%start_index(j,i)
+       surf_e = surf_def_v(0)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_def_v(0)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy2
+       ENDDO
+!
+!--    South-facing vertical default-type surfaces
+       surf_s = surf_def_v(1)%start_index(j,i)
+       surf_e = surf_def_v(1)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_def_v(1)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy2
+       ENDDO
+!
+!--    East-facing vertical default-type surfaces
+       surf_s = surf_def_v(2)%start_index(j,i)
+       surf_e = surf_def_v(2)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_def_v(2)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx2
+       ENDDO
+!
+!--    West-facing vertical default-type surfaces
+       surf_s = surf_def_v(3)%start_index(j,i)
+       surf_e = surf_def_v(3)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_def_v(3)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx2
+       ENDDO
+!
+!--    Now, for natural-type surfaces
+!--    North-facing
+       surf_s = surf_lsm_v(0)%start_index(j,i)
+       surf_e = surf_lsm_v(0)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_lsm_v(0)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy2
+       ENDDO
+!
+!--    South-facing
+       surf_s = surf_lsm_v(1)%start_index(j,i)
+       surf_e = surf_lsm_v(1)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_lsm_v(1)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy2
+       ENDDO
+!
+!--    East-facing
+       surf_s = surf_lsm_v(2)%start_index(j,i)
+       surf_e = surf_lsm_v(2)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_lsm_v(2)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx2
+       ENDDO
+!
+!--    West-facing
+       surf_s = surf_lsm_v(3)%start_index(j,i)
+       surf_e = surf_lsm_v(3)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_lsm_v(3)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx2
+       ENDDO
+!
+!--    Now, for urban-type surfaces
+!--    North-facing
+       surf_s = surf_usm_v(0)%start_index(j,i)
+       surf_e = surf_usm_v(0)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_usm_v(0)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy2
+       ENDDO
+!
+!--    South-facing
+       surf_s = surf_usm_v(1)%start_index(j,i)
+       surf_e = surf_usm_v(1)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_usm_v(1)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy2
+       ENDDO
+!
+!--    East-facing
+       surf_s = surf_usm_v(2)%start_index(j,i)
+       surf_e = surf_usm_v(2)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_usm_v(2)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx2
+       ENDDO
+!
+!--    West-facing
+       surf_s = surf_usm_v(3)%start_index(j,i)
+       surf_e = surf_usm_v(3)%end_index(j,i)
+       DO  m = surf_s, surf_e
+          k           = surf_usm_v(3)%k(m)
+          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx2
+       ENDDO
+
+
+!
+!--    Compute vertical diffusion. In case that surface fluxes have been
+!--    prescribed or computed at bottom and/or top, index k starts/ends at
+!--    nzb+2 or nzt-1, respectively. Model top is also mask if top flux 
+!--    is given. 
+       DO  k = nzb+1, nzt
+!
+!--       Determine flags to mask topography below and above. Flag 0 is 
+!--       used to mask topography in general, and flag 8 implies 
+!--       information about use_surface_fluxes. Flag 9 is used to control 
+!--       flux at model top.   
+          mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
+                               BTEST( wall_flags_0(k-1,j,i), 8 ) ) 
+          mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
+                               BTEST( wall_flags_0(k+1,j,i), 8 ) )  *          &
+                        MERGE( 1.0_wp, 0.0_wp,                                 &
+                               BTEST( wall_flags_0(k+1,j,i), 9 ) )
+          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
+                               BTEST( wall_flags_0(k,j,i), 0 ) )
+
+          tend(k,j,i) = tend(k,j,i)                                            &
+                                       + 0.5_wp * (                            &
+                                      ( kh(k,j,i) + kh(k+1,j,i) ) *            &
+                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
+                                                            * rho_air_zw(k)    &
+                                                            * mask_top         &
+                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
+                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
+                                                            * rho_air_zw(k-1)  &
+                                                            * mask_bottom      &
+                                                  ) * ddzw(k) * drho_air(k)    &
+                                                              * flag
+       ENDDO
+
+!
+!--    Vertical diffusion at horizontal walls.
+!--    TO DO: Adjust for downward facing walls and mask already in main loop
+       IF ( use_surface_fluxes )  THEN
+!
+!--       Default-type surfaces, upward-facing
+          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)
+
+             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)                    &
+                                       * ddzw(k) * drho_air(k)
+          ENDDO
+!
+!--       Default-type surfaces, downward-facing
+          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)
+
+             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)                  &
+                                       * ddzw(k) * drho_air(k)
+          ENDDO
+!
+!--       Natural-type surfaces, upward-facing
+          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) + s_flux_lsm_h_up(m)                    &
+                                       * ddzw(k) * drho_air(k)
+          ENDDO
+!
+!--       Urban-type surfaces, upward-facing
+          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) + s_flux_usm_h_up(m)                    &
+                                       * ddzw(k) * drho_air(k)
           ENDDO
        ENDIF
-
-!
-!--    Compute vertical diffusion. In case that surface fluxes have been
-!--    prescribed or computed at bottom and/or top, index k starts/ends at
-!--    nzb+2 or nzt-1, respectively.
-       DO  k = nzb_diff_s_inner(j,i), nzt_diff
-
-          tend(k,j,i) = tend(k,j,i)                                            &
-                                       + 0.5_wp * (                            &
-            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
-                                                            * rho_air_zw(k)    &
-          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                                            * rho_air_zw(k-1)  &
-                                                  ) * ddzw(k) * drho_air(k)
-       ENDDO
-
-!
-!--    Vertical diffusion at the first computational gridpoint along z-direction
-       IF ( use_surface_fluxes )  THEN
-
-          k = nzb_s_inner(j,i)+1
-
-          tend(k,j,i) = tend(k,j,i) + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )     &
-                                               * ( s(k+1,j,i)-s(k,j,i) )       &
-                                               * ddzu(k+1)                     &
-                                               * rho_air_zw(k)                 &
-                                        + s_flux_b(j,i)                        &
-                                      ) * ddzw(k) * drho_air(k)
-
-       ENDIF
-
 !
 !--    Vertical diffusion at the last computational gridpoint along z-direction
        IF ( use_top_fluxes )  THEN
-
-          k = nzt
-
-          tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                        &
-                                      - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )     &
-                                               * ( s(k,j,i)-s(k-1,j,i) )       &
-                                               * ddzu(k)                       &
-                                               * rho_air_zw(k-1)               &
-                                      ) * ddzw(k) * drho_air(k)
-
+          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)                                         &
+                           + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
+          ENDDO
        ENDIF
 
