Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 2119)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 2232)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adjustments to new topography and surface concept
 ! 
 ! Former revisions:
@@ -224,10 +224,8 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
-               momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q,    &
-               qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, &
-               sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, &
-               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
-               uswst, vsws, v, vg, vpt, vswst, w, w_subs,                      &
-               waterflux_output_conversion, zw
+               momentumflux_output_conversion, nr, p, prho, prr, pt, q,        &
+               qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s,                  &
+               sa, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, u,   &
+               ug, v, vg, vpt, w, w_subs, waterflux_output_conversion, zw
         
     USE cloud_parameters,                                                      &
@@ -236,5 +234,5 @@
     USE control_parameters,                                                    &
         ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
-                dt_3d, g, humidity, kappa, large_scale_forcing,                &
+                dt_3d, g, humidity, kappa, land_surface, large_scale_forcing,  &
                 large_scale_subsidence, max_pr_user, message_string, neutral,  &
                 microphysics_seifert, ocean, passive_scalar, simulated_time,   &
@@ -250,13 +248,10 @@
     USE indices,                                                               &
         ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
-                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
-                nzb_s_inner, nzt, nzt_diff
+                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
         
     USE kinds
     
     USE land_surface_model_mod,                                                &
-        ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
-                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
-                shf_eb, t_soil
+        ONLY:   m_soil_h, nzb_soil, nzt_soil, t_soil_h
 
     USE netcdf_interface,                                                      &
@@ -277,4 +272,7 @@
     USE statistics
 
+       USE surface_mod,                                                        &
+          ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
+
 
     IMPLICIT NONE
@@ -283,14 +281,18 @@
     INTEGER(iwp) ::  j                   !<
     INTEGER(iwp) ::  k                   !<
+    INTEGER(iwp) ::  ki                  !< 
     INTEGER(iwp) ::  k_surface_level     !<
+    INTEGER(iwp) ::  m                   !< loop variable over all horizontal wall elements 
+    INTEGER(iwp) ::  l                   !< loop variable over surface facing -- up- or downward-facing
     INTEGER(iwp) ::  nt                  !<
     INTEGER(iwp) ::  omp_get_thread_num  !<
     INTEGER(iwp) ::  sr                  !<
     INTEGER(iwp) ::  tn                  !<
-    
+
     LOGICAL ::  first  !<
     
     REAL(wp) ::  dptdz_threshold  !< 
     REAL(wp) ::  fac              !<
+    REAL(wp) ::  flag             !<
     REAL(wp) ::  height           !<
     REAL(wp) ::  pts              !<
@@ -400,15 +402,17 @@
 !--    for other horizontal averages.
        tn = 0
-
-       !$OMP PARALLEL PRIVATE( i, j, k, tn )
-!$     tn = omp_get_thread_num()
-
+       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
+       !$ tn = omp_get_thread_num()
        !$OMP DO
        DO  i = nxl, nxr
           DO  j =  nys, nyn
-             DO  k = nzb_s_inner(j,i), nzt+1
-                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
-                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
-                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
+             DO  k = nzb, nzt+1
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
+                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)  &
+                                                              * flag
+                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)  &
+                                                              * flag
+                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)  &
+                                                              * flag
              ENDDO
           ENDDO
@@ -421,7 +425,9 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                DO  k = nzb_s_inner(j,i), nzt+1
-                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
-                                      sa(k,j,i) * rmask(j,i,sr)
+                DO  k = nzb, nzt+1
+                   sums_l(k,23,tn)  = sums_l(k,23,tn) + sa(k,j,i)              &
+                                    * rmask(j,i,sr)                            &
+                                    * MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 22 ) )
                 ENDDO
              ENDDO
@@ -437,9 +443,10 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                DO  k = nzb_s_inner(j,i), nzt+1
-                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
-                                      vpt(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
-                                      q(k,j,i) * rmask(j,i,sr)
+                DO  k = nzb, nzt+1
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
+                   sums_l(k,44,tn)  = sums_l(k,44,tn) +                        &
+                                      vpt(k,j,i) * rmask(j,i,sr) * flag
+                   sums_l(k,41,tn)  = sums_l(k,41,tn) +                        &
+                                      q(k,j,i) * rmask(j,i,sr)   * flag
                 ENDDO
              ENDDO
@@ -449,10 +456,13 @@
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   DO  k = nzb_s_inner(j,i), nzt+1
-                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
-                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
-                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
+                   DO  k = nzb, nzt+1
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
+                      sums_l(k,42,tn) = sums_l(k,42,tn) +                      &
+                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) &
+                                                               * flag
+                      sums_l(k,43,tn) = sums_l(k,43,tn) + (                    &
                                       pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
-                                                          ) * rmask(j,i,sr)
+                                                          ) * rmask(j,i,sr)    &
+                                                            * flag
                    ENDDO
                 ENDDO
@@ -467,6 +477,9 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                DO  k = nzb_s_inner(j,i), nzt+1
-                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr)
+                DO  k = nzb, nzt+1
+                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i)             &
+                                    * rmask(j,i,sr)                            &
+                                    * MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 22 ) )
                 ENDDO
              ENDDO
@@ -603,39 +616,48 @@
        !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper,             &
        !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
-       !$OMP                   w2 )
-!$     tn = omp_get_thread_num()
-
+       !$OMP                   w2, flag, m, ki, l )
+       !$ tn = omp_get_thread_num()
        !$OMP DO
        DO  i = nxl, nxr
           DO  j =  nys, nyn
              sums_l_etot = 0.0_wp
-             DO  k = nzb_s_inner(j,i), nzt+1
+             DO  k = nzb, nzt+1
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
 !
 !--             Prognostic and diagnostic variables
-                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
-                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
-                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
-                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
-                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
+                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)  &
+                                                              * flag
+                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)  &
+                                                              * flag
+                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)  &
+                                                              * flag
+                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)  &
+                                                              * flag
+                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)  * flag
 
                 sums_l(k,33,tn) = sums_l(k,33,tn) + &
-                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
+                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)&
+                                                                 * flag
 
                 IF ( humidity )  THEN
                    sums_l(k,70,tn) = sums_l(k,70,tn) + &
-                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
+                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)&
+                                                                 * flag
                 ENDIF
                 IF ( passive_scalar )  THEN
                    sums_l(k,118,tn) = sums_l(k,118,tn) + &
-                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)
+                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)&
+                                                                  * flag
                 ENDIF
 !
 !--             Higher moments
 !--             (Computation of the skewness of w further below)
-                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
+                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) &
+                                                                * flag
 
                 sums_l_etot  = sums_l_etot + &
-                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + &
-                                        w(k,j,i)**2 ) * rmask(j,i,sr)
+                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 +  &
+                                        w(k,j,i)**2 )            * rmask(j,i,sr)&
+                                                                 * flag
              ENDDO
 !
@@ -647,19 +669,60 @@
 !
 !--          2D-arrays (being collected in the last column of sums_l)
-             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
-                                        us(j,i)   * rmask(j,i,sr)
-             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
-                                        usws(j,i) * rmask(j,i,sr)
-             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
-                                        vsws(j,i) * rmask(j,i,sr)
-             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
-                                        ts(j,i)   * rmask(j,i,sr)
-             IF ( humidity )  THEN
-                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
-                                            qs(j,i)   * rmask(j,i,sr)
+             IF ( surf_def_h(0)%ns >= 1 )  THEN
+                m = surf_def_h(0)%start_index(j,i)
+                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
+                                        surf_def_h(0)%us(m)   * rmask(j,i,sr)
+                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
+                                        surf_def_h(0)%usws(m) * rmask(j,i,sr)
+                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
+                                        surf_def_h(0)%vsws(m) * rmask(j,i,sr)
+                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
+                                        surf_def_h(0)%ts(m)   * rmask(j,i,sr)
+                IF ( humidity )  THEN
+                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
+                                            surf_def_h(0)%qs(m)   * rmask(j,i,sr)
+                ENDIF
+                IF ( passive_scalar )  THEN
+                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
+                                            surf_def_h(0)%ss(m)   * rmask(j,i,sr)
+                ENDIF
              ENDIF
-             IF ( passive_scalar )  THEN
-                sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
-                                            ss(j,i)   * rmask(j,i,sr)
+             IF ( surf_lsm_h%ns >= 1 )  THEN
+                m = surf_lsm_h%start_index(j,i)
+                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
+                                        surf_lsm_h%us(m)   * rmask(j,i,sr)
+                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
+                                        surf_lsm_h%usws(m) * rmask(j,i,sr)
+                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
+                                        surf_lsm_h%vsws(m) * rmask(j,i,sr)
+                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
+                                        surf_lsm_h%ts(m)   * rmask(j,i,sr)
+                IF ( humidity )  THEN
+                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
+                                            surf_lsm_h%qs(m)   * rmask(j,i,sr)
+                ENDIF
+                IF ( passive_scalar )  THEN
+                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
+                                            surf_lsm_h%ss(m)   * rmask(j,i,sr)
+                ENDIF
+             ENDIF
+             IF ( surf_usm_h%ns >= 1 )  THEN
+                m = surf_lsm_h%start_index(j,i)
+                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
+                                        surf_usm_h%us(m)   * rmask(j,i,sr)
+                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
+                                        surf_usm_h%usws(m) * rmask(j,i,sr)
+                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
+                                        surf_usm_h%vsws(m) * rmask(j,i,sr)
+                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
+                                        surf_usm_h%ts(m)   * rmask(j,i,sr)
+                IF ( humidity )  THEN
+                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
+                                            surf_usm_h%qs(m)   * rmask(j,i,sr)
+                ENDIF
+                IF ( passive_scalar )  THEN
+                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
+                                            surf_usm_h%ss(m)   * rmask(j,i,sr)
+                ENDIF
              ENDIF
           ENDDO
@@ -674,5 +737,7 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                DO  k = nzb_s_inner(j,i), nzt+1
+                DO  k = nzb, nzt+1
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
+
                    u2   = u(k,j,i)**2
                    v2   = v(k,j,i)**2
@@ -681,12 +746,16 @@
                    vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
 
-                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
-                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
-                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
+                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)    &
+                                                            * flag
+                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)    &
+                                                            * flag
+                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)    &
+                                                            * flag
 !
 !--                Perturbation energy
 
                    sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
-                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
+                                  ( ust2 + vst2 + w2 )      * rmask(j,i,sr)    &
+                                                            * flag
                 ENDDO
              ENDDO
@@ -702,5 +771,7 @@
        DO  i = nxl, nxr
           DO  j =  nys, nyn
-             DO  k = nzb_s_inner(j,i), nzt+1
+             DO  k = nzb, nzt+1
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
+
                 w2   = w(k,j,i)**2
                 ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
@@ -709,5 +780,7 @@
 
                 sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
-                                 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
+                                 + 0.5_wp * ( ust2 + vst2 + w2 )               &
+                                 * rmask(j,i,sr)                               &
+                                 * flag
              ENDDO
           ENDDO
@@ -728,6 +801,11 @@
 !--                However, this implies no error since staggered velocity 
 !--                components are zero at the walls and inside buildings.
-
-             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
+!--          Flag 23 is used to mask surface fluxes as well as model-top fluxes,
+!--          which are added further below. 
+             DO  k = nzb, nzt
+                flag = MERGE( 1.0_wp, 0.0_wp,                                  &
+                              BTEST( wall_flags_0(k,j,i), 23 ) ) *             &
+                       MERGE( 1.0_wp, 0.0_wp,                                  &
+                              BTEST( wall_flags_0(k,j,i), 9  ) )
 !
 !--             Momentum flux w"u"
@@ -739,5 +817,6 @@
                                                            ) * rmask(j,i,sr)   &
                                          * rho_air_zw(k)                       &
-                                         * momentumflux_output_conversion(k)
+                                         * momentumflux_output_conversion(k)   &
+                                         * flag
 !
 !--             Momentum flux w"v"
@@ -749,5 +828,6 @@
                                                            ) * rmask(j,i,sr)   &
                                          * rho_air_zw(k)                       &
-                                         * momentumflux_output_conversion(k)
+                                         * momentumflux_output_conversion(k)   &
+                                         * flag
 !
 !--             Heat flux w"pt"
@@ -757,5 +837,6 @@
                                                * rho_air_zw(k)                 &
                                                * heatflux_output_conversion(k) &
-                                               * ddzu(k+1) * rmask(j,i,sr)
+                                               * ddzu(k+1) * rmask(j,i,sr)     &
+                                               * flag
 
 
@@ -766,5 +847,6 @@
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( sa(k+1,j,i) - sa(k,j,i) )   &
-                                               * ddzu(k+1) * rmask(j,i,sr)
+                                               * ddzu(k+1) * rmask(j,i,sr)     &
+                                               * flag
                 ENDIF
 
@@ -777,5 +859,5 @@
                                                * rho_air_zw(k)                 &
                                                * heatflux_output_conversion(k) &
-                                               * ddzu(k+1) * rmask(j,i,sr)
+                                               * ddzu(k+1) * rmask(j,i,sr) * flag
                    sums_l(k,48,tn) = sums_l(k,48,tn)                           &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
@@ -783,5 +865,5 @@
                                                * rho_air_zw(k)                 &
                                                * waterflux_output_conversion(k)&
-                                               * ddzu(k+1) * rmask(j,i,sr)
+                                               * ddzu(k+1) * rmask(j,i,sr) * flag
 
                    IF ( cloud_physics ) THEN
@@ -792,5 +874,5 @@
                                                * rho_air_zw(k)                 &
                                                * waterflux_output_conversion(k)&
-                                               * ddzu(k+1) * rmask(j,i,sr) 
+                                               * ddzu(k+1) * rmask(j,i,sr) * flag
                    ENDIF
                 ENDIF
@@ -802,5 +884,6 @@
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                   * ( s(k+1,j,i) - s(k,j,i) )  &
-                                                  * ddzu(k+1) * rmask(j,i,sr)
+                                                  * ddzu(k+1) * rmask(j,i,sr)  &
+                                                              * flag
                 ENDIF
 
@@ -810,66 +893,179 @@
 !--          Subgridscale fluxes in the Prandtl layer
              IF ( use_surface_fluxes )  THEN
-                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
+                DO  l = 0, 1
+                   ki = MERGE( -1, 0, l == 0 )
+                   IF ( surf_def_h(l)%ns >= 1 )  THEN
+                      DO  m = surf_def_h(l)%start_index(j,i), surf_def_h(l)%end_index(j,i)
+                         k = surf_def_h(l)%k(m)
+
+                         sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + &
+                                    momentumflux_output_conversion(k+ki) * &
+                                    surf_def_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
+                         sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + &
+                                    momentumflux_output_conversion(k+ki) * &
+                                    surf_def_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
+                         sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + &
+                                    heatflux_output_conversion(k+ki) * &
+                                    surf_def_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
+                         sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + &
+                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
+                         sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
+                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
+                         IF ( ocean )  THEN
+                            sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
+                                       surf_def_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
+                         ENDIF
+                         IF ( humidity )  THEN
+                            sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) +                     &
+                                       waterflux_output_conversion(k+ki) *      &
+                                       surf_def_h(l)%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
+                            sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                   &
+                                       ( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) *     &
+                                       surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) *      &
+                                                  surf_def_h(l)%qsws(m) )                  &
+                                       * heatflux_output_conversion(k+ki)
+                            IF ( cloud_droplets )  THEN
+                               sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                &
+                                         ( 1.0_wp + 0.61_wp * q(k+ki,j,i) -     &
+                                           ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) +          &
+                                           0.61_wp * pt(k+ki,j,i) * surf_def_h(l)%qsws(m) ) &
+                                          * heatflux_output_conversion(k+ki)
+                            ENDIF
+                            IF ( cloud_physics )  THEN
+!
+!--                            Formula does not work if ql(k+ki) /= 0.0
+                               sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) +                  &
+                                          waterflux_output_conversion(k+ki) *   &
+                                          surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
+                            ENDIF
+                         ENDIF
+                         IF ( passive_scalar )  THEN
+                            sums_l(k+ki,119,tn) = sums_l(k+ki,119,tn) +                     &
+                                        surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
+                         ENDIF
+
+                      ENDDO
+
+                   ENDIF
+                ENDDO
+                IF ( surf_lsm_h%ns >= 1 )  THEN
+                   m = surf_lsm_h%start_index(j,i)
+                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
                                     momentumflux_output_conversion(nzb) * &
-                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
-                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
+                                    surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
+                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
                                     momentumflux_output_conversion(nzb) * &
-                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
-                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
+                                    surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
+                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
                                     heatflux_output_conversion(nzb) * &
-                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
-                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
+                                    surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
+                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
-                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
+                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
-                IF ( ocean )  THEN
-                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
-                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
-                ENDIF
-                IF ( humidity )  THEN
-                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
+                   IF ( ocean )  THEN
+                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
+                                       surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
+                   ENDIF
+                   IF ( humidity )  THEN
+                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
                                        waterflux_output_conversion(nzb) *      &
-                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
-                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
+                                       surf_lsm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
+                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
                                        ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
-                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
-                                                  qsws(j,i) )                  &
+                                       surf_lsm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
+                                                  surf_lsm_h%qsws(m) )                  &
                                        * heatflux_output_conversion(nzb)
-                   IF ( cloud_droplets )  THEN
-                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
+                      IF ( cloud_droplets )  THEN
+                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
-                                           ql(nzb,j,i) ) * shf(j,i) +          &
-                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) ) &
+                                           ql(nzb,j,i) ) * surf_lsm_h%shf(m) +          &
+                                           0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) ) &
                                           * heatflux_output_conversion(nzb)
+                      ENDIF
+                      IF ( cloud_physics )  THEN
+!
+!--                      Formula does not work if ql(nzb) /= 0.0
+                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
+                                          waterflux_output_conversion(nzb) *   &
+                                          surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
+                      ENDIF
                    ENDIF
-                   IF ( cloud_physics )  THEN
-!
-!--                   Formula does not work if ql(nzb) /= 0.0
-                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
+                   IF ( passive_scalar )  THEN
+                      sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
+                                        surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
+                   ENDIF
+
+
+                ENDIF
+                IF ( surf_usm_h%ns >= 1 )  THEN
+                   m = surf_usm_h%start_index(j,i)
+                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
+                                    momentumflux_output_conversion(nzb) * &
+                                    surf_usm_h%usws(m) * rmask(j,i,sr)     ! w"u"
+                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
+                                    momentumflux_output_conversion(nzb) * &
+                                    surf_usm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
+                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
+                                    heatflux_output_conversion(nzb) * &
+                                    surf_usm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
+                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
+                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
+                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
+                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
+                   IF ( ocean )  THEN
+                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
+                                       surf_usm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
+                   ENDIF
+                   IF ( humidity )  THEN
+                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
+                                       waterflux_output_conversion(nzb) *      &
+                                       surf_usm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
+                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
+                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
+                                       surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
+                                                  surf_usm_h%qsws(m) )                  &
+                                       * heatflux_output_conversion(nzb)
+                      IF ( cloud_droplets )  THEN
+                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
+                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
+                                           ql(nzb,j,i) ) * surf_usm_h%shf(m) +          &
+                                           0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) ) &
+                                          * heatflux_output_conversion(nzb)
+                      ENDIF
+                      IF ( cloud_physics )  THEN
+!
+!--                      Formula does not work if ql(nzb) /= 0.0
+                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
                                           waterflux_output_conversion(nzb) *   &
-                                          qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                                          surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
+                      ENDIF
                    ENDIF
-                ENDIF
-                IF ( passive_scalar )  THEN
-                   sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
-                                        ssws(j,i) * rmask(j,i,sr) ! w"s"
-                ENDIF
+                   IF ( passive_scalar )  THEN
+                      sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
+                                        surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
+                   ENDIF
+
+
+                ENDIF
+
              ENDIF
 
              IF ( .NOT. neutral )  THEN
-                sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
-                                    ol(j,i)  * rmask(j,i,sr) ! L
-             ENDIF
-
-
-             IF ( land_surface )  THEN
-                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + ghf_eb(j,i)
-                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + shf_eb(j,i)
-                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + qsws_eb(j,i)
-                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + qsws_liq_eb(j,i)
-                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + qsws_soil_eb(j,i)
-                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + qsws_veg_eb(j,i)
-                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn) + r_a(j,i)
-                sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ r_s(j,i)
+                IF ( surf_def_h(0)%ns >= 1 )  THEN
+                   m = surf_def_h(0)%start_index(j,i)
+                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
+                                        surf_def_h(0)%ol(m)  * rmask(j,i,sr) ! L
+                ENDIF
+                IF ( surf_lsm_h%ns >= 1 )  THEN
+                   m = surf_lsm_h%start_index(j,i)
+                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
+                                        surf_lsm_h%ol(m)  * rmask(j,i,sr) ! L
+                ENDIF
+                IF ( surf_usm_h%ns >= 1 )  THEN
+                   m = surf_usm_h%start_index(j,i)
+                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
+                                        surf_usm_h%ol(m)  * rmask(j,i,sr) ! L
+                ENDIF
              ENDIF
 
@@ -893,13 +1089,14 @@
 !--          Subgridscale fluxes at the top surface
              IF ( use_top_fluxes )  THEN
+                m = surf_def_h(2)%start_index(j,i)
                 sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
                                     momentumflux_output_conversion(nzt:nzt+1) * &
-                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
+                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
                 sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
                                     momentumflux_output_conversion(nzt:nzt+1) * &
-                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
+                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
                 sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
                                     heatflux_output_conversion(nzt:nzt+1) * &
-                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
+                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
                 sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
@@ -909,20 +1106,23 @@
                 IF ( ocean )  THEN
                    sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
-                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
+                                       surf_def_h(2)%sasws(m) * rmask(j,i,sr)  ! w"sa"
                 ENDIF
                 IF ( humidity )  THEN
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
                                        waterflux_output_conversion(nzt) *      &
-                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                                       surf_def_h(2)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
                    sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
                                        ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
-                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
-                                                             qswst(j,i) )      &
+                                       surf_def_h(2)%shf(m) +                  &
+                                       0.61_wp * pt(nzt,j,i) *    &
+                                       surf_def_h(2)%qsws(m) )      &
                                        * heatflux_output_conversion(nzt)
                    IF ( cloud_droplets )  THEN
                       sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
                                           ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
-                                            ql(nzt,j,i) ) * tswst(j,i) +       &
-                                           0.61_wp * pt(nzt,j,i) * qswst(j,i) )&
+                                            ql(nzt,j,i) ) *                    &
+                                            surf_def_h(2)%shf(m) +             &
+                                           0.61_wp * pt(nzt,j,i) *             &
+                                           surf_def_h(2)%qsws(m) )&
                                            * heatflux_output_conversion(nzt)
                    ENDIF
@@ -932,10 +1132,10 @@
                       sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
                                           waterflux_output_conversion(nzt) *   &
-                                          qswst(j,i) * rmask(j,i,sr)
+                                          surf_def_h(2)%qsws(m) * rmask(j,i,sr)
                    ENDIF
                 ENDIF
                 IF ( passive_scalar )  THEN
                    sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
-                                        sswst(j,i) * rmask(j,i,sr) ! w"s"
+                                        surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
                 ENDIF
              ENDIF
@@ -946,5 +1146,6 @@
 !--          ----  speaking the following k-loop would have to be split up and 
 !--                rearranged according to the staggered grid.
-             DO  k = nzb_s_inner(j,i), nzt
+             DO  k = nzb, nzt
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
                 ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
                                  u(k+1,j,i) - hom(k+1,1,1,sr) )
@@ -956,7 +1157,7 @@
 !--             Higher moments
                 sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
-                                                    rmask(j,i,sr)
+                                                    rmask(j,i,sr) * flag
                 sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
-                                                    rmask(j,i,sr)
+                                                    rmask(j,i,sr) * flag
 
 !
@@ -968,10 +1169,10 @@
                                        sa(k+1,j,i) - hom(k+1,1,23,sr) )
                       sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
-                                        rmask(j,i,sr)
+                                        rmask(j,i,sr) * flag
                    ENDIF
-                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *            &
-                                                       rmask(j,i,sr)
+                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *      &
+                                                       rmask(j,i,sr) * flag
                    sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
-                                                       rmask(j,i,sr)
+                                                       rmask(j,i,sr) * flag
                 ENDIF
 
@@ -985,6 +1186,7 @@
                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
                                                heatflux_output_conversion(k) * &
-                                                          rmask(j,i,sr)
-                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
+                                                          rmask(j,i,sr) * flag
+                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) &
+                                                                    * flag
 
                       IF ( .NOT. cloud_droplets )  THEN
@@ -996,14 +1198,19 @@
                          sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
                                              waterflux_output_conversion(k) *  &
-                                                             rmask(j,i,sr)
+                                                             rmask(j,i,sr)  *  &
+                                                             flag
                          sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
-                                                             rmask(j,i,sr)
+                                                             rmask(j,i,sr) *   &
+                                                             flag
                          sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
-                                                             rmask(j,i,sr)
+                                                             rmask(j,i,sr) *   &
+                                                             flag
                          IF ( microphysics_seifert )  THEN
                             sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
-                                                                rmask(j,i,sr)
+                                                                rmask(j,i,sr) *&
+                                                                flag 
                             sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
-                                                                rmask(j,i,sr)
+                                                                rmask(j,i,sr) *&
+                                                                flag
                          ENDIF
                       ENDIF
@@ -1015,5 +1222,6 @@
                          sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
                                               heatflux_output_conversion(k) *  &
-                                                             rmask(j,i,sr)
+                                                             rmask(j,i,sr)  *  &
+                                                             flag
                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
                          sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              & 
@@ -1022,5 +1230,6 @@
                                              0.61_wp * hom(k,1,4,sr) *         &
                                              sums_l(k,49,tn)                   &
-                                           ) * heatflux_output_conversion(k)
+                                           ) * heatflux_output_conversion(k) * &
+                                               flag
                       END IF
                    END IF
@@ -1033,5 +1242,5 @@
                                     s(k+1,j,i) - hom(k+1,1,117,sr) )
                    sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
-                                                       rmask(j,i,sr)
+                                                       rmask(j,i,sr) * flag
                 ENDIF
 
@@ -1042,66 +1251,132 @@
                                              ( ust**2 + vst**2 + w(k,j,i)**2 ) &
                                            * momentumflux_output_conversion(k) &
-                                             * rmask(j,i,sr)
+                                           * rmask(j,i,sr) * flag
              ENDDO
           ENDDO
        ENDDO
+       !$OMP END PARALLEL
+!
+!--    Treat land-surface quantities according to new wall model structure. 
+       IF ( land_surface )  THEN
+          tn = 0
+          !$OMP PARALLEL PRIVATE( i, j, m, tn )
+          !$ tn = omp_get_thread_num()
+          !$OMP DO
+          DO  m = 1, surf_lsm_h%ns
+             i = surf_lsm_h%i(m)
+             j = surf_lsm_h%j(m)
+        
+             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
+                  j >= nys  .AND.  j <= nyn )  THEN 
+                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf_eb(m)
+                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%shf_eb(m)
+                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h%qsws_eb(m)
+                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h%qsws_liq_eb(m)
+                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h%qsws_soil_eb(m)
+                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + surf_lsm_h%qsws_veg_eb(m)
+                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn) + surf_lsm_h%r_a(m)
+                sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ surf_lsm_h%r_s(m)
+             ENDIF
+          ENDDO
+          !$OMP END PARALLEL
+
+          tn = 0
+          !$OMP PARALLEL PRIVATE( i, j, k, m, tn )
+          !$ tn = omp_get_thread_num()
+          !$OMP DO
+          DO  m = 1, surf_lsm_h%ns
+
+             i = surf_lsm_h%i(m)            
+             j = surf_lsm_h%j(m)
+
+             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
+                  j >= nys  .AND.  j <= nyn )  THEN 
+
+                DO  k = nzb_soil, nzt_soil
+                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h%var_2d(k,m)  &
+                                      * rmask(j,i,sr)
+                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h%var_2d(k,m)  &
+                                      * rmask(j,i,sr)
+                ENDDO
+             ENDIF
+          ENDDO
+          !$OMP END PARALLEL
+       ENDIF
 !
 !--    For speed optimization fluxes which have been computed in part directly
 !--    inside the WS advection routines are treated seperatly
 !--    Momentum fluxes first:
+
+       tn = 0
+       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
+       !$ tn = omp_get_thread_num()
        IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
-         !$OMP DO
-         DO  i = nxl, nxr
-            DO  j = nys, nyn
-               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
-                  ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                &
-                                   u(k+1,j,i) - hom(k+1,1,1,sr) )
-                  vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                &
-                                   v(k+1,j,i) - hom(k+1,1,2,sr) )
-!
-!--               Momentum flux w*u*
-                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
-                                                    ( w(k,j,i-1) + w(k,j,i) )  &
-                                          * momentumflux_output_conversion(k)  &
-                                                    * ust * rmask(j,i,sr)
-!
-!--               Momentum flux w*v*
-                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
-                                                    ( w(k,j-1,i) + w(k,j,i) )  &
-                                          * momentumflux_output_conversion(k)  &
-                                                    * vst * rmask(j,i,sr)
-               ENDDO
-            ENDDO
-         ENDDO
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb, nzt
+!
+!--                Flag 23 is used to mask surface fluxes as well as model-top 
+!--                fluxes, which are added further below. 
+                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
+                                 BTEST( wall_flags_0(k,j,i), 23 ) ) *          &
+                          MERGE( 1.0_wp, 0.0_wp,                               &
+                                 BTEST( wall_flags_0(k,j,i), 9  ) )
+
+                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
+                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
+                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
+                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
+!
+!--                Momentum flux w*u*
+                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                &
+                                                     ( w(k,j,i-1) + w(k,j,i) ) &
+                                           * momentumflux_output_conversion(k) &
+                                                     * ust * rmask(j,i,sr)     &
+                                                           * flag
+!
+!--                Momentum flux w*v*
+                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                &
+                                                     ( w(k,j-1,i) + w(k,j,i) ) &
+                                           * momentumflux_output_conversion(k) &
+                                                     * vst * rmask(j,i,sr)     &
+                                                           * flag
+                ENDDO
+             ENDDO
+          ENDDO
 
        ENDIF
        IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
-         !$OMP DO
-         DO  i = nxl, nxr
-            DO  j = nys, nyn
-               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
-!
-!--               Vertical heat flux
-                  sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                 &
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb, nzt
+                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
+                                 BTEST( wall_flags_0(k,j,i), 23 ) ) *          &
+                          MERGE( 1.0_wp, 0.0_wp,                               &
+                                 BTEST( wall_flags_0(k,j,i), 9  ) )
+!
+!--                Vertical heat flux
+                   sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                &
                            ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
                              pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
                            * heatflux_output_conversion(k)                     &
-                           * w(k,j,i) * rmask(j,i,sr)
-                  IF ( humidity )  THEN
-                     pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +            &
+                           * w(k,j,i) * rmask(j,i,sr) * flag
+                   IF ( humidity )  THEN
+                      pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +           &
                                       q(k+1,j,i) - hom(k+1,1,41,sr) )
-                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
+                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *     &
                                        waterflux_output_conversion(k) *        &
-                                       rmask(j,i,sr)
-                  ENDIF
-                  IF ( passive_scalar )  THEN
-                     pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +            &
+                                       rmask(j,i,sr) * flag
+                   ENDIF
+                   IF ( passive_scalar )  THEN
+                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +           &
                                       s(k+1,j,i) - hom(k+1,1,117,sr) )
-                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *     &
-                                        rmask(j,i,sr)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDDO
+                      sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *    &
+                                        rmask(j,i,sr) * flag
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
 
        ENDIF
@@ -1126,5 +1401,6 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
+                DO  k = nzb+1, nzt
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
                    sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
@@ -1133,8 +1409,8 @@
                 + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
                             - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
-                + w(k,j,i)**2                                        )
+                + w(k,j,i)**2                                        ) * flag
 
                    sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
-                                               * ( p(k,j,i) + p(k+1,j,i) )
+                                               * ( p(k,j,i) + p(k+1,j,i) ) * flag
 
                 ENDDO
@@ -1164,14 +1440,17 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
+                DO  k = nzb+1, nzt
+
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
                    sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
                  - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
-                                                                ) * ddzw(k)
+                                                                ) * ddzw(k)    &
+                                                                  * flag
 
                    sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
-                                                                )
+                                                                )  * flag
 
                 ENDDO
@@ -1191,5 +1470,6 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
+                DO  k = nzb+1, nzt
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 !
 !--                Subgrid horizontal heat fluxes u"pt", v"pt"
@@ -1199,5 +1479,5 @@
                                                * rho_air_zw(k)                 &
                                                * heatflux_output_conversion(k) &
-                                                 * ddx * rmask(j,i,sr)
+                                                 * ddx * rmask(j,i,sr) * flag
                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
@@ -1205,5 +1485,5 @@
                                                * rho_air_zw(k)                 &
                                                * heatflux_output_conversion(k) &
-                                                 * ddy * rmask(j,i,sr)
+                                                 * ddy * rmask(j,i,sr) * flag
 !
 !--                Resolved horizontal heat fluxes u*pt*, v*pt*
@@ -1212,5 +1492,6 @@
                                     * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
                                                  pt(k,j,i)   - hom(k,1,4,sr) ) &
-                                               * heatflux_output_conversion(k)
+                                               * heatflux_output_conversion(k) &
+                                               * flag
                    pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
                                     pt(k,j,i)   - hom(k,1,4,sr) )
@@ -1219,5 +1500,6 @@
                                     * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
                                                  pt(k,j,i)   - hom(k,1,4,sr) ) &
-                                               * heatflux_output_conversion(k)
+                                               * heatflux_output_conversion(k) &
+                                               * flag
                 ENDDO
              ENDDO
@@ -1279,41 +1561,30 @@
        ENDIF
 
+       tn = 0
        !$OMP PARALLEL PRIVATE( i, j, k, tn )
-!$     tn = omp_get_thread_num()
-       IF ( land_surface )  THEN
-          !$OMP DO
-          DO  i = nxl, nxr
-             DO  j =  nys, nyn
-                DO  k = nzb_soil, nzt_soil
-                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
-                                      * rmask(j,i,sr)
-                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
-                                      * rmask(j,i,sr)
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
+       !$ tn = omp_get_thread_num()       
        IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
           !$OMP DO
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt+1
+                DO  k = nzb+1, nzt+1
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
                    sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
-                                       * rmask(j,i,sr)
+                                       * rmask(j,i,sr) * flag
                    sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
-                                       * rmask(j,i,sr)
+                                       * rmask(j,i,sr) * flag
                    sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
-                                       * rmask(j,i,sr)
+                                       * rmask(j,i,sr) * flag
                    sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
-                                       * rmask(j,i,sr)
+                                       * rmask(j,i,sr) * flag
                    sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
-                                       * rmask(j,i,sr)
+                                       * rmask(j,i,sr) * flag
                    sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
-                                       * rmask(j,i,sr)
+                                       * rmask(j,i,sr) * flag
                    sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
-                                       * rmask(j,i,sr)
+                                       * rmask(j,i,sr) * flag
                    sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
-                                       * rmask(j,i,sr)
+                                       * rmask(j,i,sr) * flag
                 ENDDO
              ENDDO
