source: palm/trunk/SOURCE/flow_statistics.f90 @ 18

 SUBROUTINE flow_statistics

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: flow_statistics.f90 4 2007-02-13 11:33:16Z raasch $
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.41  2006/08/04 14:37:50  raasch
! Error removed in non-parallel part (sums_l)
!
! Revision 1.1  1997/08/11 06:15:17  raasch
! Initial revision
!
!
! Description:
! ------------
! Compute average profiles and further average flow quantities for the different
! 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.
!------------------------------------------------------------------------------!

    USE arrays_3d
    USE cloud_parameters
    USE cpulog
    USE grid_variables
    USE indices
    USE interfaces
    USE pegrid
    USE statistics
    USE control_parameters

    IMPLICIT NONE

    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
    LOGICAL ::  first
    REAL    ::  height, pts, sums_l_eper, sums_l_etot, ust, ust2, u2, vst, &
                vst2, v2, w2, z_i(2)
    REAL    ::  sums_ll(nzb:nzt+1,2)


    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )

!
!-- To be on the safe side, check whether flow_statistics has already been
!-- called once after the current time step
    IF ( flow_statistics_called )  THEN
       IF ( myid == 0 )  PRINT*, '+++ WARNING: flow_statistics is called two', &
                                 ' times within one timestep'
       CALL local_stop
    ENDIF

!
!-- Compute statistics for each (sub-)region
    DO  sr = 0, statistic_regions

!
!--    Initialize (local) summation array
       sums_l = 0.0

!
!--    Store sums that have been computed in other subroutines in summation
!--    array
       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
!--    WARNING: next line still has to be adjusted for OpenMP 
       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
       sums_l(nzb+9,var_sum,0)  = sums_divold_l(sr)  ! old divergence from pres
       sums_l(nzb+10,var_sum,0) = sums_divnew_l(sr)  ! new divergence from pres
!--    WARNING: next four lines still may have to be adjusted for OpenMP 
       sums_l(nzb:nzb+2,var_sum-1,0)    = sums_up_fraction_l(1,1:3,sr)! upstream
       sums_l(nzb+3:nzb+5,var_sum-1,0)  = sums_up_fraction_l(2,1:3,sr)! parts
       sums_l(nzb+6:nzb+8,var_sum-1,0)  = sums_up_fraction_l(3,1:3,sr)! from
       sums_l(nzb+9:nzb+11,var_sum-1,0) = sums_up_fraction_l(4,1:3,sr)! spline

!
!--    Horizontally averaged profiles of horizontal velocities and temperature.
!--    They must have been computed before, because they are already required
!--    for other horizontal averages.
       tn = 0
       !$OMP PARALLEL PRIVATE( i, j, k, tn )
#if defined( __lcmuk )
       tn = omp_get_thread_num()
#else
!$     tn = omp_get_thread_num()
#endif

       !$OMP DO
       DO  i = nxl, nxr
          DO  j =  nys, nyn
             DO  k = nzb_s_outer(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)
             ENDDO
          ENDDO
       ENDDO

!
!--    Horizontally averaged profiles of virtual potential temperature,
!--    total water content, specific humidity and liquid water potential
!--    temperature
       IF ( moisture )  THEN
          !$OMP DO
          DO  i = nxl, nxr
             DO  j =  nys, nyn
                DO  k = nzb_s_outer(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)
                ENDDO
             ENDDO
          ENDDO
          IF ( cloud_physics )  THEN
             !$OMP DO
             DO  i = nxl, nxr
                DO  j =  nys, nyn
                   DO  k = nzb_s_outer(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) + ( &
                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
                                                          ) * rmask(j,i,sr)
                   ENDDO
                ENDDO
             ENDDO
          ENDIF
       ENDIF

!
!--    Horizontally averaged profiles of passive scalar
       IF ( passive_scalar )  THEN
          !$OMP DO
          DO  i = nxl, nxr
             DO  j =  nys, nyn
                DO  k = nzb_s_outer(j,i), nzt+1
                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
                ENDDO
             ENDDO
          ENDDO
       ENDIF
       !$OMP END PARALLEL

!
!--    Summation of thread sums
       IF ( threads_per_task > 1 )  THEN
          DO  i = 1, threads_per_task-1
             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
             IF ( moisture )  THEN
                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
                IF ( cloud_physics )  THEN
                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
                ENDIF
             ENDIF
             IF ( passive_scalar )  THEN
                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
             ENDIF
          ENDDO
       ENDIF

#if defined( __parallel )
!
!--    Compute total sum from local sums
       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
                           MPI_SUM, comm2d, ierr )
       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
                           MPI_SUM, comm2d, ierr )
       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
                           MPI_SUM, comm2d, ierr )
       IF ( moisture ) THEN
          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
                              MPI_REAL, MPI_SUM, comm2d, ierr )
          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
                              MPI_REAL, MPI_SUM, comm2d, ierr )
          IF ( cloud_physics ) THEN
             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
                                 MPI_REAL, MPI_SUM, comm2d, ierr )
             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
                                 MPI_REAL, MPI_SUM, comm2d, ierr )
          ENDIF
       ENDIF

       IF ( passive_scalar )  THEN
          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
                              MPI_REAL, MPI_SUM, comm2d, ierr )
       ENDIF
#else
       sums(:,1) = sums_l(:,1,0)
       sums(:,2) = sums_l(:,2,0)
       sums(:,4) = sums_l(:,4,0)
       IF ( moisture ) THEN
          sums(:,44) = sums_l(:,44,0)
          sums(:,41) = sums_l(:,41,0)
          IF ( cloud_physics ) THEN
             sums(:,42) = sums_l(:,42,0)
             sums(:,43) = sums_l(:,43,0)
          ENDIF
       ENDIF
       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
#endif

!
!--    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)
       hom(:,1,1,sr) = sums(:,1)             ! u
       hom(:,1,2,sr) = sums(:,2)             ! v
       hom(:,1,4,sr) = sums(:,4)             ! pt

!
!--    Humidity and cloud parameters
       IF ( moisture ) THEN
          sums(:,44) = sums(:,44) / ngp_2dh_outer(:,sr)
          sums(:,41) = sums(:,41) / ngp_2dh_outer(:,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)
             hom(:,1,42,sr) = sums(:,42)             ! qv
             hom(:,1,43,sr) = sums(:,43)             ! pt
          ENDIF
       ENDIF

!
!--    Passive scalar
       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) / ngp_2dh_outer(:,sr)

!
!--    Horizontally averaged profiles of the remaining prognostic variables,
!--    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 
!--    ----  speaking the following k-loop would have to be split up and 
!--          rearranged according to the staggered grid.
       tn = 0
#if defined( __lcmuk )
       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
       tn = omp_get_thread_num()
#else
       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
!$     tn = omp_get_thread_num()
#endif
       !$OMP DO
       DO  i = nxl, nxr
          DO  j =  nys, nyn
             sums_l_etot = 0.0
             sums_l_eper = 0.0
             DO  k = nzb_s_outer(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
                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)

!
!--             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)
!
!--             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)
             ENDDO
!
!--          Total and perturbation energy for the total domain (being
!--          collected in the last column of sums_l). Summation of these
!--          quantities is seperated from the previous loop in order to
!--          allow vectorization of that loop.
             sums_l(nzb+4,var_sum,tn) = sums_l(nzb+4,var_sum,tn) + sums_l_etot
             sums_l(nzb+5,var_sum,tn) = sums_l(nzb+5,var_sum,tn) + sums_l_eper
!
!--          2D-arrays (being collected in the last column of sums_l)
             sums_l(nzb,var_sum,tn)   = sums_l(nzb,var_sum,tn) +   &
                                        us(j,i)   * rmask(j,i,sr)
             sums_l(nzb+1,var_sum,tn) = sums_l(nzb+1,var_sum,tn) + &
                                        usws(j,i) * rmask(j,i,sr)
             sums_l(nzb+2,var_sum,tn) = sums_l(nzb+2,var_sum,tn) + &
                                        vsws(j,i) * rmask(j,i,sr)
             sums_l(nzb+3,var_sum,tn) = sums_l(nzb+3,var_sum,tn) + &
                                        ts(j,i)   * rmask(j,i,sr)
          ENDDO
       ENDDO

!
!--    Horizontally averaged profiles of the vertical fluxes
       !$OMP DO
       DO  i = nxl, nxr
          DO  j = nys, nyn
!
!--          Subgridscale fluxes (without Prandtl layer from k=nzb, 
!--          oterwise from k=nzb+1)
!--          NOTE: for simplicity, nzb_diff_s_outer 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
!
!--             Momentum flux w"u"
                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * (                   &
                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
                                                           ) * (               &
                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
                                                               ) * rmask(j,i,sr)
!
!--             Momentum flux w"v"
                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * (                   &
                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
                                                           ) * (               &
                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
                                                               ) * rmask(j,i,sr)
!
!--             Heat flux w"pt"
                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
                                               * ddzu(k+1) * rmask(j,i,sr)


!
!--             Buoyancy flux, water flux (humidity flux) w"q"
                IF ( moisture ) THEN
                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
                                               * ddzu(k+1) * rmask(j,i,sr)
                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
                                               * ( q(k+1,j,i) - q(k,j,i) )     &
                                               * ddzu(k+1) * rmask(j,i,sr)
                   IF ( cloud_physics ) THEN
                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
                                               * ddzu(k+1) * rmask(j,i,sr) 
                   ENDIF
                ENDIF

!
!--             Passive scalar flux
                IF ( passive_scalar )  THEN
                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
                                               * ( q(k+1,j,i) - q(k,j,i) )     &
                                               * ddzu(k+1) * rmask(j,i,sr)
                ENDIF

             ENDDO

!
!--          Subgridscale fluxes in the Prandtl layer
             IF ( use_surface_fluxes )  THEN
                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
                                    0.0 * rmask(j,i,sr)           ! u"pt"
                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
                                    0.0 * rmask(j,i,sr)           ! v"pt"
                IF ( moisture )  THEN
                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
                   IF ( cloud_physics )  THEN
                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
                                          ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
                                          shf(j,i) + 0.61 * pt(nzb,j,i) * &
                                                     qsws(j,i)            &
                                                              )
!
!--                   Formula does not work if ql(nzb) /= 0.0
                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
                                          qsws(j,i) * rmask(j,i,sr)
                   ENDIF
                ENDIF
                IF ( passive_scalar )  THEN
                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
                ENDIF
             ENDIF

!
!--          Resolved fluxes (can be computed for all horizontal points)
!--          NOTE: for simplicity, nzb_s_outer 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
                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) )
                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 * &
                                                    rmask(j,i,sr)
                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
                                                    rmask(j,i,sr)

!
!--             Buoyancy flux, water flux, humidity flux and liquid water
!--             content
                IF ( moisture )  THEN
                   pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
                                 vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                   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 *                                            &
                             ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) &
                             + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
                      sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
                                                          rmask(j,i,sr)
                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
                                                          rmask(j,i,sr)
                   ENDIF
                ENDIF

!
!--             Passive scalar flux
                IF ( passive_scalar )  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

!
!--             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)
          
             ENDDO
          ENDDO
       ENDDO

!
!--    Divergence of vertical flux of resolved scale energy and pressure
!--    fluctuations. First calculate the products, then the divergence.
!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
       IF ( hom(nzb+1,2,55,0) /= 0.0 )  THEN

          sums_ll = 0.0  ! local array

          !$OMP DO
          DO  i = nxl, nxr
             DO  j = nys, nyn
                DO  k = nzb_s_outer(j,i)+1, nzt

                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
                           - 2.0 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
                           ) )**2                                          &
                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
                           - 2.0 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
                           ) )**2                                          &
                   + w(k,j,i)**2                                  )

                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
                                               * ( p(k,j,i) + p(k+1,j,i) )

                ENDDO
             ENDDO
          ENDDO
          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
          sums_ll(nzt+1,1) = 0.0
          sums_ll(0,2)     = 0.0
          sums_ll(nzt+1,2) = 0.0

          DO  k = nzb_s_outer(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)
          ENDDO
          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)

       ENDIF

!
!--    Divergence of vertical flux of SGS TKE
       IF ( hom(nzb+1,2,57,0) /= 0.0 )  THEN

          !$OMP DO
          DO  i = nxl, nxr
             DO  j = nys, nyn
                DO  k = nzb_s_outer(j,i)+1, nzt

                   sums_l(k,57,tn) = sums_l(k,57,tn) + (                       &
                   (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)

                ENDDO
             ENDDO
          ENDDO
          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)

       ENDIF

!
!--    Horizontal heat fluxes (subgrid, resolved, total). 
!--    Do it only, if profiles shall be plotted. 
       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN

          !$OMP DO
          DO  i = nxl, nxr
             DO  j = nys, nyn
                DO  k = nzb_s_outer(j,i)+1, nzt
!
!--                Subgrid horizontal heat fluxes u"pt", v"pt"
                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
                                                 * ddx * rmask(j,i,sr)
                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
                                                 * ddy * rmask(j,i,sr)
!
!--                Resolved horizontal heat fluxes u*pt*, v*pt*
                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
                                                 pt(k,j,i)   - hom(k,1,4,sr) )
                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
                                 pt(k,j,i)   - hom(k,1,4,sr) )
                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
                                                 pt(k,j,i)   - hom(k,1,4,sr) )
                ENDDO
             ENDDO
          ENDDO
!
!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
          sums(nzb,58) = 0.0
          sums(nzb,59) = 0.0
          sums(nzb,60) = 0.0
          sums(nzb,61) = 0.0
          sums(nzb,62) = 0.0
          sums(nzb,63) = 0.0

       ENDIF
       !$OMP END PARALLEL

!
!--    Summation of thread sums
       IF ( threads_per_task > 1 )  THEN
          DO  i = 1, threads_per_task-1
             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
             sums_l(:,45:var_sum,0) = sums_l(:,45:var_sum,0) + &
                                      sums_l(:,45:var_sum,i)
          ENDDO
       ENDIF

#if defined( __parallel )
!
!--    Compute total sum from local sums
       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
                           MPI_SUM, comm2d, ierr )
#else
       sums = sums_l(:,:,0)
#endif

!
!--    Final values are obtained by division by the total number of grid points 
!--    used for summation. After that store profiles.
!--    Profiles:
       DO  k = nzb, nzt+1
          sums(k,:var_sum-2)      = sums(k,:var_sum-2) / ngp_2dh_outer(k,sr)
       ENDDO
!--    Upstream-parts
       sums(nzb:nzb+11,var_sum-1) = sums(nzb:nzb+11,var_sum-1) / ngp_3d(sr)
!--    u* and so on
!--    As sums(nzb:nzb+3,var_sum) are full 2D arrays (us, usws, vsws, ts) whose
!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
!--    above the topography, they are being divided by ngp_2dh(sr)
       sums(nzb:nzb+3,var_sum)    = sums(nzb:nzb+3,var_sum)    / &
                                    ngp_2dh(sr)
!--    eges, e*
       sums(nzb+4:nzb+5,var_sum)  = sums(nzb+4:nzb+5,var_sum)  / &
                                    ngp_3d_inner(sr)
!--    Old and new divergence
       sums(nzb+9:nzb+10,var_sum) = sums(nzb+9:nzb+10,var_sum) / &
                                    ngp_3d_inner(sr)

!
!--    Collect horizontal average in hom.
!--    Compute deduced averages (e.g. total heat flux)
       hom(:,1,3,sr)  = sums(:,3)      ! w
       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
       hom(:,1,9,sr)  = sums(:,9)      ! km
       hom(:,1,10,sr) = sums(:,10)     ! kh
       hom(:,1,11,sr) = sums(:,11)     ! l
       hom(:,1,12,sr) = sums(:,12)     ! w"u"
       hom(:,1,13,sr) = sums(:,13)     ! w*u*
       hom(:,1,14,sr) = sums(:,14)     ! w"v"
       hom(:,1,15,sr) = sums(:,15)     ! w*v*
       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
                                       ! profiles 23-29 left empty for initial 
                                       ! profiles
       hom(:,1,30,sr) = sums(:,30)     ! u*2
       hom(:,1,31,sr) = sums(:,31)     ! v*2
       hom(:,1,32,sr) = sums(:,32)     ! w*2
       hom(:,1,33,sr) = sums(:,33)     ! pt*2
       hom(:,1,34,sr) = sums(:,34)     ! e*
       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
       hom(:,1,37,sr) = sums(:,37)     ! w*e*
       hom(:,1,38,sr) = sums(:,38)     ! w*3
       hom(:,1,39,sr) = sums(:,38) / ( sums(:,32) + 1E-20 )**1.5    ! Sw
       hom(:,1,40,sr) = sums(:,40)     ! p
       hom(:,1,45,sr) = sums(:,45)     ! w"q"
       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
       hom(:,1,54,sr) = sums(:,54)     ! ql
       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
       hom(:,1,57,sr) = sums(:,57)     ! w"e/dz
       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t

       hom(:,1,var_hom-1,sr) = sums(:,var_sum-1)
                                       ! upstream-parts u_x, u_y, u_z, v_x,
                                       ! v_y, usw. (in last but one profile)
       hom(:,1,var_hom,sr) = sums(:,var_sum)  
                                       ! u*, w'u', w'v', t* (in last profile)

!
!--    Determine the boundary layer height using two different schemes.
!--    First scheme: Starting from the Earth's surface, look for the first 
!--    relative minimum of the total heat flux. The corresponding height is
!--    accepted as the boundary layer height, if it is less than 1.5 times the
!--    height where the heat flux becomes negative for the first time.
!--    NOTE: This criterion is still capable of improving!
       z_i(1) = 0.0
       first = .TRUE.
       DO  k = nzb, nzt-1
          IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
             first = .FALSE.
             height = zw(k)
          ENDIF
          IF ( hom(k,1,18,sr) < 0.0  .AND. &
               hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
             IF ( zw(k) < 1.5 * height )  THEN
                z_i(1) = zw(k)
             ELSE
                z_i(1) = height
             ENDIF
             EXIT
          ENDIF
       ENDDO

!
!--    Second scheme: Starting from the top model boundary, look for the first 
!--    characteristic kink in the temperature profile, where the originally 
!--    stable stratification notably weakens.
       z_i(2) = 0.0
       DO  k = nzt-1, nzb+1, -1
          IF ( ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) > &
               2.0 * ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) )  THEN
             z_i(2) = zu(k)
             EXIT
          ENDIF
       ENDDO

       hom(nzb+6,1,var_hom,sr) = z_i(1)
       hom(nzb+7,1,var_hom,sr) = z_i(2)

!
!--    Computation of both the characteristic vertical velocity and
!--    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
          hom(nzb+8,1,var_hom,sr)  = ( g / hom(nzb+1,1,4,sr) * &
                                       hom(nzb,1,18,sr) * z_i(1) )**0.333333333
!--       so far this only works if Prandtl layer is used
          hom(nzb+11,1,var_hom,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,var_hom,sr)
       ELSE
          hom(nzb+8,1,var_hom,sr)  = 0.0
          hom(nzb+11,1,var_hom,sr) = 0.0
       ENDIF

    ENDDO    ! loop of the subregions

!
!-- Calculate additional statistics provided by the user interface
    CALL user_statistics

!
!-- If required, sum up horizontal averages for subsequent time averaging
    IF ( do_sum )  THEN
       IF ( average_count_pr == 0 )  hom_sum = 0.0
       hom_sum = hom_sum + hom(:,1,:,:)
       average_count_pr = average_count_pr + 1
       do_sum = .FALSE.
    ENDIF

!
!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
!-- This flag is reset after each time step in time_integration.
    flow_statistics_called = .TRUE.

    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )


 END SUBROUTINE flow_statistics