Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 110)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 132)
@@ -4,4 +4,7 @@
 ! Actual revisions:
 ! -----------------
+! Vertical profiles now based on nzb_s_inner; they are divided by 
+! ngp_2dh_s_inner (scalars, procucts of scalars and velocity components) and 
+! ngp_2dh (staggered velocity components and their products), respectively.
 ! 
 !
@@ -49,10 +52,9 @@
 ! user-defined (sub-)regions. The region indexed 0 is the total model domain.
 !
-! NOTE: For simplicity, nzb_s_outer and nzb_diff_s_outer are being used as a
-! ----  lower vertical index for k-loops for all variables so that regardless
-! of the variable and its respective staggered grid always the same number of
-! grid points is used for 2D averages. The disadvantage: depending on the
-! variable, up to one grid layer adjacent to the (vertical walls of the)
-! topography is missed out by this simplification.
+! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
+! ----  lower vertical index for k-loops for all variables, although strictly 
+! speaking the k-loops would have to be split up according to the staggered 
+! grid. However, this implies no error since staggered velocity components are 
+! zero at the walls and inside buildings.
 !------------------------------------------------------------------------------!
 
@@ -124,5 +126,5 @@
        DO  i = nxl, nxr
           DO  j =  nys, nyn
-             DO  k = nzb_s_outer(j,i), nzt+1
+             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)
@@ -138,5 +140,5 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                DO  k = nzb_s_outer(j,i), nzt+1
+                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)
@@ -154,5 +156,5 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                DO  k = nzb_s_outer(j,i), nzt+1
+                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)
@@ -166,5 +168,5 @@
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   DO  k = nzb_s_outer(j,i), nzt+1
+                   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)
@@ -184,5 +186,5 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                DO  k = nzb_s_outer(j,i), nzt+1
+                DO  k = nzb_s_inner(j,i), nzt+1
                    sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
                 ENDDO
@@ -265,7 +267,7 @@
 !--    Final values are obtained by division by the total number of grid points 
 !--    used for summation. After that store profiles.
-       sums(:,1) = sums(:,1) / ngp_2dh_outer(:,sr)
-       sums(:,2) = sums(:,2) / ngp_2dh_outer(:,sr)
-       sums(:,4) = sums(:,4) / ngp_2dh_outer(:,sr)
+       sums(:,1) = sums(:,1) / ngp_2dh(sr)
+       sums(:,2) = sums(:,2) / ngp_2dh(sr)
+       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
        hom(:,1,1,sr) = sums(:,1)             ! u
        hom(:,1,2,sr) = sums(:,2)             ! v
@@ -275,5 +277,5 @@
 !--    Salinity
        IF ( ocean )  THEN
-          sums(:,23) = sums(:,23) / ngp_2dh_outer(:,sr)
+          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
           hom(:,1,23,sr) = sums(:,23)             ! sa
        ENDIF
@@ -282,11 +284,11 @@
 !--    Humidity and cloud parameters
        IF ( humidity ) THEN
-          sums(:,44) = sums(:,44) / ngp_2dh_outer(:,sr)
-          sums(:,41) = sums(:,41) / ngp_2dh_outer(:,sr)
+          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
+          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
           hom(:,1,44,sr) = sums(:,44)             ! vpt
           hom(:,1,41,sr) = sums(:,41)             ! qv (q)
           IF ( cloud_physics ) THEN
-             sums(:,42) = sums(:,42) / ngp_2dh_outer(:,sr)
-             sums(:,43) = sums(:,43) / ngp_2dh_outer(:,sr)
+             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
+             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
              hom(:,1,42,sr) = sums(:,42)             ! qv
              hom(:,1,43,sr) = sums(:,43)             ! pt
@@ -296,5 +298,6 @@
 !
 !--    Passive scalar
-       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) / ngp_2dh_outer(:,sr)
+       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /  &
+            ngp_2dh_s_inner(:,sr)                    ! s (q)
 
 !
@@ -302,7 +305,9 @@
 !--    variances, the total and the perturbation energy (single values in last
 !--    column of sums_l) and some diagnostic quantities.
-!--    NOTE: for simplicity, nzb_s_outer is used below, although strictly 
+!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly 
 !--    ----  speaking the following k-loop would have to be split up and 
 !--          rearranged according to the staggered grid.
+!--          However, this implies no error since staggered velocity components
+!--          are zero at the walls and inside buildings.
        tn = 0
 #if defined( __intel_openmp_bug )
@@ -319,5 +324,5 @@
              sums_l_etot = 0.0
              sums_l_eper = 0.0
-             DO  k = nzb_s_outer(j,i), nzt+1
+             DO  k = nzb_s_inner(j,i), nzt+1
                 u2   = u(k,j,i)**2
                 v2   = v(k,j,i)**2
@@ -382,8 +387,11 @@
 !--          Subgridscale fluxes (without Prandtl layer from k=nzb, 
 !--          oterwise from k=nzb+1)
-!--          NOTE: for simplicity, nzb_diff_s_outer is used below, although
+!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
 !--          ----  strictly speaking the following k-loop would have to be
 !--                split up according to the staggered grid.
-             DO  k = nzb_diff_s_outer(j,i)-1, nzt_diff
+!--                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
 !
 !--             Momentum flux w"u"
@@ -528,8 +536,8 @@
 !
 !--          Resolved fluxes (can be computed for all horizontal points)
-!--          NOTE: for simplicity, nzb_s_outer is used below, although strictly 
+!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly 
 !--          ----  speaking the following k-loop would have to be split up and 
 !--                rearranged according to the staggered grid.
-             DO  k = nzb_s_outer(j,i), nzt
+             DO  k = nzb_s_inner(j,i), nzt
                 ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
                               u(k+1,j,i) - hom(k+1,1,1,sr) )
@@ -634,5 +642,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_outer(j,i)+1, nzt
+                DO  k = nzb_s_inner(j,i)+1, nzt
 
                    sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
@@ -656,5 +664,5 @@
           sums_ll(nzt+1,2) = 0.0
 
-          DO  k = nzb_s_outer(j,i)+1, nzt
+          DO  k = nzb_s_inner(j,i)+1, nzt
              sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
              sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
@@ -674,5 +682,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_outer(j,i)+1, nzt
+                DO  k = nzb_s_inner(j,i)+1, nzt
 
                    sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
@@ -701,5 +709,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_outer(j,i)+1, nzt
+                DO  k = nzb_s_inner(j,i)+1, nzt
 !
 !--                Subgrid horizontal heat fluxes u"pt", v"pt"
@@ -773,5 +781,18 @@
 !--    Profiles:
        DO  k = nzb, nzt+1
-          sums(k,:pr_palm-2)      = sums(k,:pr_palm-2) / ngp_2dh_outer(k,sr)
+          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
+          sums(k,9:11)            = sums(k,9:11)        / ngp_2dh_s_inner(k,sr)
+          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
+          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
+          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
+          sums(k,33)              = sums(k,33)          / ngp_2dh_s_inner(k,sr)
+          sums(k,34:39)           = sums(k,34:39)       / ngp_2dh(sr)
+          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
+          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
+          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
+          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
+          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
+          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
+          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
 !--    Upstream-parts
@@ -785,5 +806,5 @@
 !--    eges, e*
        sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
-                                    ngp_3d_inner(sr)
+                                    ngp_3d(sr)
 !--    Old and new divergence
        sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
@@ -795,5 +816,5 @@
              sums(k,pr_palm+1:pr_palm+max_pr_user) = &
                                     sums(k,pr_palm+1:pr_palm+max_pr_user) / &
-                                    ngp_2dh_outer(k,sr)
+                                    ngp_2dh_s_inner(k,sr)
           ENDDO
        ENDIF
