Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 625)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 667)
@@ -4,5 +4,14 @@
 ! Current revisions:
 ! -----------------
-! 
+! When advection is computed with ws-scheme, turbulent fluxes are already
+! computed in the respective advection routines and buffered in arrays
+! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the
+! numerics and to avoid unphysical kinks near the surface. 
+! So some if requests has to be done to dicern between fluxes from ws-scheme
+! other advection schemes. 
+! Furthermore the computation of z_i is only done if the heat flux exceeds a
+! minimum value. This affects only simulations of a neutral boundary layer and
+! is due to reasons of computations in the advection scheme.
+!
 !
 ! Former revisions:
@@ -102,5 +111,4 @@
     REAL    ::  sums_ll(nzb:nzt+1,2)
 
-
     CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
 
@@ -133,4 +141,41 @@
        sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
 
+!
+!--    Copy the turbulent quantities, evaluated in the advection routines to 
+!--    the local array sums_l() for further computations
+       IF ( ws_scheme_mom )  THEN
+!       
+!--       Boundary condition for u'u' and v'v', because below the surface no
+!--       computation for these quantities is done.
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                sums_us2_ws_l(nzb_u_inner(j,i),sr) =                          & 
+                    sums_us2_ws_l(nzb_u_inner(j,i)+1,sr)
+                sums_vs2_ws_l(nzb_v_inner(j,i),sr) =                          &  
+                    sums_vs2_ws_l(nzb_v_inner(j,i)+1,sr)
+             ENDDO
+          ENDDO
+!          
+!--       Swap the turbulent quantities evaluated in advec_ws.
+          sums_l(:,13,0) = sums_wsus_ws_l(:,sr)       ! w*u* 
+          sums_l(:,15,0) = sums_wsvs_ws_l(:,sr)       ! w*v* 
+          sums_l(:,30,0) = sums_us2_ws_l(:,sr)        ! u*2 
+          sums_l(:,31,0) = sums_vs2_ws_l(:,sr)        ! v*2 
+          sums_l(:,32,0) = sums_ws2_ws_l(:,sr)        ! w*2 
+          sums_l(:,34,0) = sums_l(:,34,0) + 0.5 *                             & 
+                (sums_us2_ws_l(:,sr) + sums_vs2_ws_l(:,sr)                    &
+                + sums_ws2_ws_l(:,sr))                      ! e*
+          DO  k = nzb, nzt
+             sums_l(nzb+5,pr_palm,0) = sums_l(nzb+5,pr_palm,0) + 0.5 * (      &
+                sums_us2_ws_l(k,sr) + sums_vs2_ws_l(k,sr) +                   &
+                sums_ws2_ws_l(k,sr))
+          ENDDO
+       ENDIF
+       IF ( ws_scheme_sca )  THEN
+          sums_l(:,17,0) = sums_wspts_ws_l(:,sr)      ! w*pts* from advec_s_ws
+          IF ( ocean ) sums_l(:,66,0) = sums_wssas_ws_l(:,sr) ! w*sa*
+          IF ( humidity  .OR.  passive_scalar ) sums_l(:,49,0) =              &
+                                                   sums_wsqs_ws_l(:,sr) !w*q* 
+       ENDIF
 ! 
 !--    Horizontally averaged profiles of horizontal velocities and temperature.
@@ -138,4 +183,5 @@
 !--    for other horizontal averages.
        tn = 0
+
        !$OMP PARALLEL PRIVATE( i, j, k, tn )
 #if defined( __intel_openmp_bug )
@@ -215,5 +261,4 @@
        ENDIF
        !$OMP END PARALLEL
-
 !
 !--    Summation of thread sums
@@ -304,4 +349,5 @@
        hom(:,1,2,sr) = sums(:,2)             ! v
        hom(:,1,4,sr) = sums(:,4)             ! pt
+
 
 !
@@ -354,11 +400,5 @@
           DO  j =  nys, nyn
              sums_l_etot = 0.0
-             sums_l_eper = 0.0
              DO  k = nzb_s_inner(j,i), nzt+1
-                u2   = u(k,j,i)**2
-                v2   = v(k,j,i)**2
-                w2   = w(k,j,i)**2
-                ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
-                vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
 !
 !--             Prognostic and diagnostic variables
@@ -369,9 +409,4 @@
                 sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
 
-!
-!--             Variances
-                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,33,tn) = sums_l(k,33,tn) + &
                                   ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
@@ -381,17 +416,8 @@
                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
                 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) * w2 * &
-                                                    rmask(j,i,sr)
-!
-!--             Perturbation energy
-                sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * ( ust2 + vst2 + w2 ) &
-                                                    * rmask(j,i,sr)
+
                 sums_l_etot  = sums_l_etot + &
-                                        0.5 * ( u2 + v2 + w2 ) * rmask(j,i,sr)
-                sums_l_eper  = sums_l_eper + &
-                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
+                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
+                                        w(k,j,i)**2 ) * rmask(j,i,sr)
              ENDDO
 !
@@ -401,5 +427,4 @@
 !--          allow vectorization of that loop.
              sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
-             sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) + sums_l_eper
 !
 !--          2D-arrays (being collected in the last column of sums_l)
@@ -420,5 +445,58 @@
 
 !
+!--    Computation of statistics when ws-scheme is not used. Else these 
+!--    quantities are evaluated in the advection routines.
+       IF ( .NOT. ws_scheme_mom )  THEN
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                sums_l_eper = 0.0
+                DO  k = nzb_s_inner(j,i), nzt+1
+                   u2   = u(k,j,i)**2
+                   v2   = v(k,j,i)**2
+                   w2   = w(k,j,i)**2
+                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
+                   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)
+!
+!--   Higher moments
+!--  (Computation of the skewness of w further below)
+                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
+                                                    rmask(j,i,sr)
+!
+!--             Perturbation energy
+
+                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
+                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
+                   sums_l_eper  = sums_l_eper + &
+                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
+
+                ENDDO
+                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
+                     + sums_l_eper
+             ENDDO
+          ENDDO
+       ELSE
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                DO  k = nzb_s_inner(j,i), nzt + 1
+                   w2   = w(k,j,i)**2
+!
+!--                Higher moments
+!--                (Computation of the skewness of w further below)
+                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
+                                                    rmask(j,i,sr)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+
+!
 !--    Horizontally averaged profiles of the vertical fluxes
+
        !$OMP DO
        DO  i = nxl, nxr
@@ -587,24 +665,5 @@
                 pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
                               pt(k+1,j,i) - hom(k+1,1,4,sr) )
-!
-!--             Momentum flux w*u*
-                sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                     &
-                                                    ( w(k,j,i-1) + w(k,j,i) ) &
-                                                    * ust * rmask(j,i,sr)
-!
-!--             Momentum flux w*v*
-                sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                     &
-                                                    ( w(k,j-1,i) + w(k,j,i) ) &
-                                                    * vst * rmask(j,i,sr)
-!
-!--             Heat flux w*pt*
-!--             The following formula (comment line, not executed) does not
-!--             work if applied to subregions
-!                sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 *                     &
-!                                                    ( pt(k,j,i)+pt(k+1,j,i) ) &
-!                                                    * w(k,j,i) * rmask(j,i,sr)
-                sums_l(k,17,tn) = sums_l(k,17,tn) + pts * w(k,j,i) * &
-                                                    rmask(j,i,sr)
-!
+
 !--             Higher moments
                 sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
@@ -617,8 +676,10 @@
 !--             but so far there is no other suitable place to calculate)
                 IF ( ocean )  THEN
-                   pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
+                   IF( .NOT. ws_scheme_sca )  THEN
+                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
                                  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) * &
+                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
                                                        rmask(j,i,sr)
+                   ENDIF
                    sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
                                                        rmask(j,i,sr)
@@ -635,8 +696,5 @@
                    sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
                                                        rmask(j,i,sr)
-                   pts = 0.5 * ( 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) * &
-                                                       rmask(j,i,sr)
+
                    IF ( cloud_physics  .OR.  cloud_droplets )  THEN
                       pts = 0.5 *                                            &
@@ -652,5 +710,5 @@
 !
 !--             Passive scalar flux
-                IF ( passive_scalar )  THEN
+                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca ))  THEN
                    pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
                                  q(k+1,j,i) - hom(k+1,1,41,sr) )
@@ -661,11 +719,62 @@
 !
 !--             Energy flux w*e*
-                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *           &
-                                              ( ust**2 + vst**2 + w(k,j,i)**2 )&
-                                              * rmask(j,i,sr)
-          
+!--             has to be adjusted
+                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
+                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
+                                             * rmask(j,i,sr)
              ENDDO
           ENDDO
        ENDDO
+!-     for reasons of speed optimization the loop is splitted, to avoid if-else
+!-     statements inside the loop
+!-     Fluxes which have been computed in part directly inside the advection routines
+!-     treated seperatly.
+!-     First treat the momentum fluxes
+       IF ( .NOT. ws_scheme_mom )  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 * ( u(k,j,i)   - hom(k,1,1,sr) + &
+                              u(k+1,j,i) - hom(k+1,1,1,sr) )
+                  vst = 0.5 * ( 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 *                   &
+                                                    ( w(k,j,i-1) + w(k,j,i) ) &
+                                                    * ust * rmask(j,i,sr)
+!
+!--               Momentum flux w*v*
+                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
+                                                    ( w(k,j-1,i) + w(k,j,i) ) &
+                                                    * vst * rmask(j,i,sr)
+               ENDDO
+            ENDDO
+         ENDDO
+
+       ENDIF
+       IF ( .NOT. ws_scheme_sca )  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 * &
+                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
+                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
+                           * w(k,j,i) * rmask(j,i,sr)
+                  IF ( humidity )  THEN
+                     pts = 0.5 * ( 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) * &
+                                                      rmask(j,i,sr)
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDDO
+
+       ENDIF
+
 
 !
@@ -814,4 +923,5 @@
 
 #if defined( __parallel )
+
 !
 !--    Compute total sum from local sums
@@ -843,4 +953,5 @@
           sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
+
 !--    Upstream-parts
        sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
@@ -868,5 +979,4 @@
           ENDDO
        ENDIF
-
 !
 !--    Collect horizontal average in hom.
@@ -934,5 +1044,5 @@
                                        ! upstream-parts u_x, u_y, u_z, v_x,
                                        ! v_y, usw. (in last but one profile)
-       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)  
+       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
 
@@ -951,11 +1061,14 @@
        z_i(1) = 0.0
        first = .TRUE.
+
        IF ( ocean )  THEN
           DO  k = nzt, nzb+1, -1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
+             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
              IF ( hom(k,1,18,sr) < 0.0  .AND. &
+                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
                   hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
                 IF ( zw(k) < 1.5 * height )  THEN
@@ -969,9 +1082,11 @@
        ELSE
           DO  k = nzb, nzt-1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
+             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0  .AND. &
+             IF ( hom(k,1,18,sr) < 0.0  .AND. & 
+                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
                   hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
                 IF ( zw(k) < 1.5 * height )  THEN
@@ -1026,5 +1141,6 @@
 !--    the characteristic convective boundary layer temperature.
 !--    The horizontal average at nzb+1 is input for the average temperature.
-       IF ( hom(nzb,1,18,sr) > 0.0  .AND.  z_i(1) /= 0.0 )  THEN
+       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
+           .AND.  z_i(1) /= 0.0 )  THEN
           hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
                                        hom(nzb,1,18,sr) *      &
