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

Last change on this file since 1267 was 1258, checked in by raasch, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 120.5 KB
RevLine 
[1221]1#if ! defined( __openacc )
[1]2 SUBROUTINE flow_statistics
3
[1036]4!--------------------------------------------------------------------------------!
5! This file is part of PALM.
6!
7! PALM is free software: you can redistribute it and/or modify it under the terms
8! of the GNU General Public License as published by the Free Software Foundation,
9! either version 3 of the License, or (at your option) any later version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
18! Copyright 1997-2012  Leibniz University Hannover
19!--------------------------------------------------------------------------------!
20!
[254]21! Current revisions:
[1]22! -----------------
[1008]23!
[1258]24!
[1008]25! Former revisions:
26! -----------------
27! $Id: flow_statistics.f90 1258 2013-11-08 16:09:09Z heinze $
28!
[1258]29! 1257 2013-11-08 15:18:40Z raasch
30! openacc "end parallel" replaced by "end parallel loop"
31!
[1242]32! 1241 2013-10-30 11:36:58Z heinze
33! Output of ug and vg
34!
[1222]35! 1221 2013-09-10 08:59:13Z raasch
36! ported for openACC in separate #else branch
37!
[1182]38! 1179 2013-06-14 05:57:58Z raasch
39! comment for profile 77 added
40!
[1116]41! 1115 2013-03-26 18:16:16Z hoffmann
42! ql is calculated by calc_liquid_water_content
43!
[1112]44! 1111 2013-03-08 23:54:10Z raasch
45! openACC directive added
46!
[1054]47! 1053 2012-11-13 17:11:03Z hoffmann
[1112]48! additions for two-moment cloud physics scheme:
[1054]49! +nr, qr, qc, prr
50!
[1037]51! 1036 2012-10-22 13:43:42Z raasch
52! code put under GPL (PALM 3.9)
53!
[1008]54! 1007 2012-09-19 14:30:36Z franke
[1007]55! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
56! turbulent fluxes of WS-scheme
57! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
58! droplets at nzb and nzt added
[700]59!
[802]60! 801 2012-01-10 17:30:36Z suehring
61! Calculation of turbulent fluxes in advec_ws is now thread-safe.
62!
[744]63! 743 2011-08-18 16:10:16Z suehring
64! Calculation of turbulent fluxes with WS-scheme only for the whole model
65! domain, not for user-defined subregions.
66!
[710]67! 709 2011-03-30 09:31:40Z raasch
68! formatting adjustments
69!
[700]70! 699 2011-03-22 17:52:22Z suehring
[699]71! Bugfix in calculation of vertical velocity skewness. The added absolute value
72! avoid negative values in the root. Negative values of w'w' can occur at the
73! top or bottom of the model domain due to degrading the order of advection
74! scheme. Furthermore the calculation will be the same for all advection
75! schemes.
[1221]76!, tend
[697]77! 696 2011-03-18 07:03:49Z raasch
78! Bugfix: Summation of Wicker-Skamarock scheme fluxes and variances for all
79! threads
80!
[679]81! 678 2011-02-02 14:31:56Z raasch
82! Bugfix in calculation of the divergence of vertical flux of resolved scale
83! energy, pressure fluctuations, and flux of pressure fluctuation itself
84!
[674]85! 673 2011-01-18 16:19:48Z suehring
86! Top bc for the horizontal velocity variances added for ocean runs.
87! Setting the corresponding bottom bc moved to advec_ws.
88!
[668]89! 667 2010-12-23 12:06:00Z suehring/gryschka
90! When advection is computed with ws-scheme, turbulent fluxes are already
91! computed in the respective advection routines and buffered in arrays
92! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the
93! numerics and to avoid unphysical kinks near the surface.
94! So some if requests has to be done to dicern between fluxes from ws-scheme
95! other advection schemes.
96! Furthermore the computation of z_i is only done if the heat flux exceeds a
97! minimum value. This affects only simulations of a neutral boundary layer and
98! is due to reasons of computations in the advection scheme.
99!
[625]100! 624 2010-12-10 11:46:30Z heinze
101! Calculation of q*2 added
102!
[623]103! 622 2010-12-10 08:08:13Z raasch
104! optional barriers included in order to speed up collective operations
105!
[392]106! 388 2009-09-23 09:40:33Z raasch
[388]107! Vertical profiles of potential density and hydrostatic pressure are
108! calculated.
[343]109! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
110! end.
[291]111! Temperature gradient criterion for estimating the boundary layer height
112! replaced by the gradient criterion of Sullivan et al. (1998).
[254]113! Output of messages replaced by message handling routine.
[1]114!
[198]115! 197 2008-09-16 15:29:03Z raasch
116! Spline timeseries splptx etc. removed, timeseries w'u', w'v', w'q' (k=0)
117! added,
118! bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr)
119! (like other scalars)
120!
[139]121! 133 2007-11-20 10:10:53Z letzel
122! Vertical profiles now based on nzb_s_inner; they are divided by
123! ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered
124! velocity components and their products, procucts of scalars and velocity
125! components), respectively.
126!
[110]127! 106 2007-08-16 14:30:26Z raasch
128! Prescribed momentum fluxes at the top surface are used,
129! profiles for w*p* and w"e are calculated
130!
[98]131! 97 2007-06-21 08:23:15Z raasch
132! Statistics for ocean version (salinity, density) added,
133! calculation of z_i and Deardorff velocity scale adjusted to be used with
134! the ocean version
135!
[90]136! 87 2007-05-22 15:46:47Z raasch
137! Two more arguments added to user_statistics, which is now also called for
138! user-defined profiles,
139! var_hom and var_sum renamed pr_palm
140!
[83]141! 82 2007-04-16 15:40:52Z raasch
142! Cpp-directive lcmuk changed to intel_openmp_bug
143!
[77]144! 75 2007-03-22 09:54:05Z raasch
145! Collection of time series quantities moved from routine flow_statistics to
146! here, routine user_statistics is called for each statistic region,
147! moisture renamed humidity
148!
[39]149! 19 2007-02-23 04:53:48Z raasch
[77]150! fluxes at top modified (tswst, qswst)
[39]151!
[3]152! RCS Log replace by Id keyword, revision history cleaned up
153!
[1]154! Revision 1.41  2006/08/04 14:37:50  raasch
155! Error removed in non-parallel part (sums_l)
156!
157! Revision 1.1  1997/08/11 06:15:17  raasch
158! Initial revision
159!
160!
161! Description:
162! ------------
163! Compute average profiles and further average flow quantities for the different
164! user-defined (sub-)regions. The region indexed 0 is the total model domain.
165!
[132]166! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
167! ----  lower vertical index for k-loops for all variables, although strictly
168! speaking the k-loops would have to be split up according to the staggered
169! grid. However, this implies no error since staggered velocity components are
170! zero at the walls and inside buildings.
[1]171!------------------------------------------------------------------------------!
172
173    USE arrays_3d
174    USE cloud_parameters
[709]175    USE control_parameters
[1]176    USE cpulog
177    USE grid_variables
178    USE indices
179    USE interfaces
180    USE pegrid
181    USE statistics
182
183    IMPLICIT NONE
184
185    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
186    LOGICAL ::  first
[291]187    REAL    ::  dptdz_threshold, height, pts, sums_l_eper, sums_l_etot, ust, &
188                ust2, u2, vst, vst2, v2, w2, z_i(2)
189    REAL    ::  dptdz(nzb+1:nzt+1)
[1]190    REAL    ::  sums_ll(nzb:nzt+1,2)
191
192    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
193
[1221]194    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, usws, v, vsws, w )
195
[1]196!
197!-- To be on the safe side, check whether flow_statistics has already been
198!-- called once after the current time step
199    IF ( flow_statistics_called )  THEN
[254]200
[274]201       message_string = 'flow_statistics is called two times within one ' // &
202                        'timestep'
[254]203       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
[1007]204
[1]205    ENDIF
206
207!
208!-- Compute statistics for each (sub-)region
209    DO  sr = 0, statistic_regions
210
211!
212!--    Initialize (local) summation array
213       sums_l = 0.0
214
215!
216!--    Store sums that have been computed in other subroutines in summation
217!--    array
218       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
219!--    WARNING: next line still has to be adjusted for OpenMP
220       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
[87]221       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
222       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
[1]223
[667]224!
225!--    Copy the turbulent quantities, evaluated in the advection routines to
226!--    the local array sums_l() for further computations
[743]227       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
[696]228
[1007]229!
[673]230!--       According to the Neumann bc for the horizontal velocity components,
231!--       the corresponding fluxes has to satisfiy the same bc.
232          IF ( ocean )  THEN
[801]233             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
[1007]234             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
[673]235          ENDIF
[696]236
237          DO  i = 0, threads_per_task-1
[1007]238!
[696]239!--          Swap the turbulent quantities evaluated in advec_ws.
[801]240             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
241             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
242             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
243             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
244             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
[696]245             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        & 
[801]246                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
247                                sums_ws2_ws_l(:,i) )    ! e*
[696]248             DO  k = nzb, nzt
[801]249                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
250                                                      sums_us2_ws_l(k,i) +  &
251                                                      sums_vs2_ws_l(k,i) +  &
252                                                      sums_ws2_ws_l(k,i) )
[696]253             ENDDO
[667]254          ENDDO
[696]255
[667]256       ENDIF
[696]257
[743]258       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
[696]259
260          DO  i = 0, threads_per_task-1
[801]261             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
262             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
[696]263             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
[801]264                                                   sums_wsqs_ws_l(:,i) !w*q*
[696]265          ENDDO
266
[667]267       ENDIF
[305]268!
[1]269!--    Horizontally averaged profiles of horizontal velocities and temperature.
270!--    They must have been computed before, because they are already required
271!--    for other horizontal averages.
272       tn = 0
[667]273
[1]274       !$OMP PARALLEL PRIVATE( i, j, k, tn )
[82]275#if defined( __intel_openmp_bug )
[1]276       tn = omp_get_thread_num()
277#else
278!$     tn = omp_get_thread_num()
279#endif
280
281       !$OMP DO
282       DO  i = nxl, nxr
283          DO  j =  nys, nyn
[132]284             DO  k = nzb_s_inner(j,i), nzt+1
[1]285                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
286                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
287                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
288             ENDDO
289          ENDDO
290       ENDDO
291
292!
[96]293!--    Horizontally averaged profile of salinity
294       IF ( ocean )  THEN
295          !$OMP DO
296          DO  i = nxl, nxr
297             DO  j =  nys, nyn
[132]298                DO  k = nzb_s_inner(j,i), nzt+1
[96]299                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
300                                      sa(k,j,i) * rmask(j,i,sr)
301                ENDDO
302             ENDDO
303          ENDDO
304       ENDIF
305
306!
[1]307!--    Horizontally averaged profiles of virtual potential temperature,
308!--    total water content, specific humidity and liquid water potential
309!--    temperature
[75]310       IF ( humidity )  THEN
[1]311          !$OMP DO
312          DO  i = nxl, nxr
313             DO  j =  nys, nyn
[132]314                DO  k = nzb_s_inner(j,i), nzt+1
[1]315                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
316                                      vpt(k,j,i) * rmask(j,i,sr)
317                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
318                                      q(k,j,i) * rmask(j,i,sr)
319                ENDDO
320             ENDDO
321          ENDDO
322          IF ( cloud_physics )  THEN
323             !$OMP DO
324             DO  i = nxl, nxr
325                DO  j =  nys, nyn
[132]326                   DO  k = nzb_s_inner(j,i), nzt+1
[1]327                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
328                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
329                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
330                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
331                                                          ) * rmask(j,i,sr)
332                   ENDDO
333                ENDDO
334             ENDDO
335          ENDIF
336       ENDIF
337
338!
339!--    Horizontally averaged profiles of passive scalar
340       IF ( passive_scalar )  THEN
341          !$OMP DO
342          DO  i = nxl, nxr
343             DO  j =  nys, nyn
[132]344                DO  k = nzb_s_inner(j,i), nzt+1
[1]345                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
346                ENDDO
347             ENDDO
348          ENDDO
349       ENDIF
350       !$OMP END PARALLEL
351!
352!--    Summation of thread sums
353       IF ( threads_per_task > 1 )  THEN
354          DO  i = 1, threads_per_task-1
355             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
356             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
357             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
[96]358             IF ( ocean )  THEN
359                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
360             ENDIF
[75]361             IF ( humidity )  THEN
[1]362                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
363                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
364                IF ( cloud_physics )  THEN
365                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
366                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
367                ENDIF
368             ENDIF
369             IF ( passive_scalar )  THEN
370                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
371             ENDIF
372          ENDDO
373       ENDIF
374
375#if defined( __parallel )
376!
377!--    Compute total sum from local sums
[622]378       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]379       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
380                           MPI_SUM, comm2d, ierr )
[622]381       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]382       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
383                           MPI_SUM, comm2d, ierr )
[622]384       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]385       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
386                           MPI_SUM, comm2d, ierr )
[96]387       IF ( ocean )  THEN
[622]388          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[96]389          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
390                              MPI_REAL, MPI_SUM, comm2d, ierr )
391       ENDIF
[75]392       IF ( humidity ) THEN
[622]393          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]394          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
395                              MPI_REAL, MPI_SUM, comm2d, ierr )
[622]396          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]397          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
398                              MPI_REAL, MPI_SUM, comm2d, ierr )
399          IF ( cloud_physics ) THEN
[622]400             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]401             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
402                                 MPI_REAL, MPI_SUM, comm2d, ierr )
[622]403             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]404             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
405                                 MPI_REAL, MPI_SUM, comm2d, ierr )
406          ENDIF
407       ENDIF
408
409       IF ( passive_scalar )  THEN
[622]410          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]411          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
412                              MPI_REAL, MPI_SUM, comm2d, ierr )
413       ENDIF
414#else
415       sums(:,1) = sums_l(:,1,0)
416       sums(:,2) = sums_l(:,2,0)
417       sums(:,4) = sums_l(:,4,0)
[96]418       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
[75]419       IF ( humidity ) THEN
[1]420          sums(:,44) = sums_l(:,44,0)
421          sums(:,41) = sums_l(:,41,0)
422          IF ( cloud_physics ) THEN
423             sums(:,42) = sums_l(:,42,0)
424             sums(:,43) = sums_l(:,43,0)
425          ENDIF
426       ENDIF
427       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
428#endif
429
430!
431!--    Final values are obtained by division by the total number of grid points
432!--    used for summation. After that store profiles.
[132]433       sums(:,1) = sums(:,1) / ngp_2dh(sr)
434       sums(:,2) = sums(:,2) / ngp_2dh(sr)
435       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
[1]436       hom(:,1,1,sr) = sums(:,1)             ! u
437       hom(:,1,2,sr) = sums(:,2)             ! v
438       hom(:,1,4,sr) = sums(:,4)             ! pt
439
[667]440
[1]441!
[96]442!--    Salinity
443       IF ( ocean )  THEN
[132]444          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
[96]445          hom(:,1,23,sr) = sums(:,23)             ! sa
446       ENDIF
447
448!
[1]449!--    Humidity and cloud parameters
[75]450       IF ( humidity ) THEN
[132]451          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
452          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
[1]453          hom(:,1,44,sr) = sums(:,44)             ! vpt
454          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
455          IF ( cloud_physics ) THEN
[132]456             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
457             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
[1]458             hom(:,1,42,sr) = sums(:,42)             ! qv
459             hom(:,1,43,sr) = sums(:,43)             ! pt
460          ENDIF
461       ENDIF
462
463!
464!--    Passive scalar
[132]465       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /  &
466            ngp_2dh_s_inner(:,sr)                    ! s (q)
[1]467
468!
469!--    Horizontally averaged profiles of the remaining prognostic variables,
470!--    variances, the total and the perturbation energy (single values in last
471!--    column of sums_l) and some diagnostic quantities.
[132]472!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
[1]473!--    ----  speaking the following k-loop would have to be split up and
474!--          rearranged according to the staggered grid.
[132]475!--          However, this implies no error since staggered velocity components
476!--          are zero at the walls and inside buildings.
[1]477       tn = 0
[82]478#if defined( __intel_openmp_bug )
[1]479       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
480       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
481       tn = omp_get_thread_num()
482#else
483       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
484!$     tn = omp_get_thread_num()
485#endif
486       !$OMP DO
487       DO  i = nxl, nxr
488          DO  j =  nys, nyn
489             sums_l_etot = 0.0
[132]490             DO  k = nzb_s_inner(j,i), nzt+1
[1]491!
492!--             Prognostic and diagnostic variables
493                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
494                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
495                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
496                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
497                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
498
499                sums_l(k,33,tn) = sums_l(k,33,tn) + &
500                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
[624]501
502                IF ( humidity )  THEN
503                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
504                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
505                ENDIF
[1007]506
[699]507!
508!--             Higher moments
509!--             (Computation of the skewness of w further below)
510                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
[667]511
[1]512                sums_l_etot  = sums_l_etot + &
[667]513                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
514                                        w(k,j,i)**2 ) * rmask(j,i,sr)
[1]515             ENDDO
516!
517!--          Total and perturbation energy for the total domain (being
518!--          collected in the last column of sums_l). Summation of these
519!--          quantities is seperated from the previous loop in order to
520!--          allow vectorization of that loop.
[87]521             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
[1]522!
523!--          2D-arrays (being collected in the last column of sums_l)
[87]524             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +   &
[1]525                                        us(j,i)   * rmask(j,i,sr)
[87]526             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
[1]527                                        usws(j,i) * rmask(j,i,sr)
[87]528             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
[1]529                                        vsws(j,i) * rmask(j,i,sr)
[87]530             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
[1]531                                        ts(j,i)   * rmask(j,i,sr)
[197]532             IF ( humidity )  THEN
533                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
534                                            qs(j,i)   * rmask(j,i,sr)
535             ENDIF
[1]536          ENDDO
537       ENDDO
538
539!
[667]540!--    Computation of statistics when ws-scheme is not used. Else these
541!--    quantities are evaluated in the advection routines.
[743]542       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 )  THEN
[667]543          !$OMP DO
544          DO  i = nxl, nxr
545             DO  j =  nys, nyn
546                sums_l_eper = 0.0
547                DO  k = nzb_s_inner(j,i), nzt+1
548                   u2   = u(k,j,i)**2
549                   v2   = v(k,j,i)**2
550                   w2   = w(k,j,i)**2
551                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
552                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
553
554                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
555                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
556                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
557!
558!--             Perturbation energy
559
560                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
561                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
562                   sums_l_eper  = sums_l_eper + &
563                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
564
565                ENDDO
566                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
567                     + sums_l_eper
568             ENDDO
569          ENDDO
570       ENDIF
[1241]571
[667]572!
[1]573!--    Horizontally averaged profiles of the vertical fluxes
[667]574
[1]575       !$OMP DO
576       DO  i = nxl, nxr
577          DO  j = nys, nyn
578!
579!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
580!--          oterwise from k=nzb+1)
[132]581!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
[1]582!--          ----  strictly speaking the following k-loop would have to be
583!--                split up according to the staggered grid.
[132]584!--                However, this implies no error since staggered velocity
585!--                components are zero at the walls and inside buildings.
586
587             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
[1]588!
589!--             Momentum flux w"u"
590                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * (                   &
591                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
592                                                           ) * (               &
593                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
594                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
595                                                               ) * rmask(j,i,sr)
596!
597!--             Momentum flux w"v"
598                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * (                   &
599                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
600                                                           ) * (               &
601                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
602                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
603                                                               ) * rmask(j,i,sr)
604!
605!--             Heat flux w"pt"
606                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
607                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
608                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
609                                               * ddzu(k+1) * rmask(j,i,sr)
610
611
612!
[96]613!--             Salinity flux w"sa"
614                IF ( ocean )  THEN
615                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
616                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
617                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
618                                               * ddzu(k+1) * rmask(j,i,sr)
619                ENDIF
620
621!
[1]622!--             Buoyancy flux, water flux (humidity flux) w"q"
[75]623                IF ( humidity ) THEN
[1]624                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
625                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
626                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
627                                               * ddzu(k+1) * rmask(j,i,sr)
628                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
629                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
630                                               * ( q(k+1,j,i) - q(k,j,i) )     &
631                                               * ddzu(k+1) * rmask(j,i,sr)
[1007]632
[1]633                   IF ( cloud_physics ) THEN
634                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
635                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
636                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
637                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
638                                               * ddzu(k+1) * rmask(j,i,sr) 
639                   ENDIF
640                ENDIF
641
642!
643!--             Passive scalar flux
644                IF ( passive_scalar )  THEN
645                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
646                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
647                                               * ( q(k+1,j,i) - q(k,j,i) )     &
648                                               * ddzu(k+1) * rmask(j,i,sr)
649                ENDIF
650
651             ENDDO
652
653!
654!--          Subgridscale fluxes in the Prandtl layer
655             IF ( use_surface_fluxes )  THEN
656                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
657                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
658                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
659                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
660                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
661                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
662                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
663                                    0.0 * rmask(j,i,sr)           ! u"pt"
664                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
665                                    0.0 * rmask(j,i,sr)           ! v"pt"
[96]666                IF ( ocean )  THEN
667                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
668                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
669                ENDIF
[75]670                IF ( humidity )  THEN
[1]671                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
672                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
[1007]673                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
674                                       ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
675                                       shf(j,i) + 0.61 * pt(nzb,j,i) * &
676                                                  qsws(j,i) )
677                   IF ( cloud_droplets )  THEN
678                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (        &
679                                         ( 1.0 + 0.61 * q(nzb,j,i) -   &
680                                           ql(nzb,j,i) ) * shf(j,i) +  &
681                                           0.61 * pt(nzb,j,i) * qsws(j,i) )
682                   ENDIF
[1]683                   IF ( cloud_physics )  THEN
684!
685!--                   Formula does not work if ql(nzb) /= 0.0
686                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
687                                          qsws(j,i) * rmask(j,i,sr)
688                   ENDIF
689                ENDIF
690                IF ( passive_scalar )  THEN
691                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
692                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
693                ENDIF
694             ENDIF
695
696!
[19]697!--          Subgridscale fluxes at the top surface
698             IF ( use_top_fluxes )  THEN
[550]699                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
[102]700                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
[550]701                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
[102]702                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
[550]703                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
[19]704                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
[550]705                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
[19]706                                    0.0 * rmask(j,i,sr)           ! u"pt"
[550]707                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
708                                    0.0 * rmask(j,i,sr)           ! v"pt"
709
[96]710                IF ( ocean )  THEN
711                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
712                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
713                ENDIF
[75]714                IF ( humidity )  THEN
[19]715                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
[388]716                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
[1007]717                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (             &
718                                       ( 1.0 + 0.61 * q(nzt,j,i) ) *     &
719                                       tswst(j,i) + 0.61 * pt(nzt,j,i) * &
720                                                  qswst(j,i) )
721                   IF ( cloud_droplets )  THEN
722                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
723                                          ( 1.0 + 0.61 * q(nzt,j,i) -     &
724                                            ql(nzt,j,i) ) * tswst(j,i) +  &
725                                            0.61 * pt(nzt,j,i) * qswst(j,i) )
726                   ENDIF
[19]727                   IF ( cloud_physics )  THEN
728!
729!--                   Formula does not work if ql(nzb) /= 0.0
730                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
731                                          qswst(j,i) * rmask(j,i,sr)
732                   ENDIF
733                ENDIF
734                IF ( passive_scalar )  THEN
735                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
[388]736                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
[19]737                ENDIF
738             ENDIF
739
740!
[1]741!--          Resolved fluxes (can be computed for all horizontal points)
[132]742!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
[1]743!--          ----  speaking the following k-loop would have to be split up and
744!--                rearranged according to the staggered grid.
[132]745             DO  k = nzb_s_inner(j,i), nzt
[1]746                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
747                              u(k+1,j,i) - hom(k+1,1,1,sr) )
748                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
749                              v(k+1,j,i) - hom(k+1,1,2,sr) )
750                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
751                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
[667]752
[1]753!--             Higher moments
754                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
755                                                    rmask(j,i,sr)
756                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
757                                                    rmask(j,i,sr)
758
759!
[96]760!--             Salinity flux and density (density does not belong to here,
[97]761!--             but so far there is no other suitable place to calculate)
[96]762                IF ( ocean )  THEN
[743]763                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
[667]764                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
[96]765                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
[667]766                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
[96]767                                                       rmask(j,i,sr)
[667]768                   ENDIF
[96]769                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
770                                                       rmask(j,i,sr)
[388]771                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
772                                                       rmask(j,i,sr)
[96]773                ENDIF
774
775!
[1053]776!--             Buoyancy flux, water flux, humidity flux, liquid water
777!--             content, rain drop concentration and rain water content
[75]778                IF ( humidity )  THEN
[1007]779                   IF ( cloud_physics .OR. cloud_droplets )  THEN
[1053]780                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +        &
[1007]781                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
782                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
[1]783                                                          rmask(j,i,sr)
[1053]784                      IF ( .NOT. cloud_droplets )  THEN
[1115]785                         pts = 0.5 *                                           &
786                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
787                              hom(k,1,42,sr) +                                 &
788                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
[1053]789                              hom(k+1,1,42,sr) )
[1115]790                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
[1053]791                                                             rmask(j,i,sr)
792                         IF ( icloud_scheme == 0  )  THEN
[1115]793                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
[1053]794                                                                rmask(j,i,sr)
[1115]795                            sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *    &
[1053]796                                                                rmask(j,i,sr)
[1115]797                            IF ( precipitation )  THEN
798                               sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * &
799                                                                   rmask(j,i,sr)
800                               sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * &
801                                                                   rmask(j,i,sr)
802                               sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *&
803                                                                   rmask(j,i,sr)
804                            ENDIF
[1053]805                         ELSE
[1115]806                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
[1053]807                                                                rmask(j,i,sr)
808                         ENDIF
809                      ELSE
[1115]810                         sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *       &
[1053]811                                                             rmask(j,i,sr)
812                      ENDIF
[1007]813                   ELSE
814                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
815                         pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
816                                       vpt(k+1,j,i) - hom(k+1,1,44,sr) )
817                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
818                                                             rmask(j,i,sr)
819                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
[1053]820                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) *  &
[1007]821                                             sums_l(k,17,tn) +      &
822                                          0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
823                      END IF
824                   END IF
[1]825                ENDIF
826!
827!--             Passive scalar flux
[743]828                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca   &
829                     .OR. sr /= 0 ) )  THEN
[1]830                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
831                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
832                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
833                                                       rmask(j,i,sr)
834                ENDIF
835
836!
837!--             Energy flux w*e*
[667]838!--             has to be adjusted
839                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
840                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
841                                             * rmask(j,i,sr)
[1]842             ENDDO
843          ENDDO
844       ENDDO
[709]845!
846!--    For speed optimization fluxes which have been computed in part directly
847!--    inside the WS advection routines are treated seperatly
848!--    Momentum fluxes first:
[743]849       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
[667]850         !$OMP DO
851         DO  i = nxl, nxr
852            DO  j = nys, nyn
853               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
854                  ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
855                              u(k+1,j,i) - hom(k+1,1,1,sr) )
856                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
857                              v(k+1,j,i) - hom(k+1,1,2,sr) )
[1007]858!
[667]859!--               Momentum flux w*u*
860                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
861                                                    ( w(k,j,i-1) + w(k,j,i) ) &
862                                                    * ust * rmask(j,i,sr)
863!
864!--               Momentum flux w*v*
865                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
866                                                    ( w(k,j-1,i) + w(k,j,i) ) &
867                                                    * vst * rmask(j,i,sr)
868               ENDDO
869            ENDDO
870         ENDDO
[1]871
[667]872       ENDIF
[743]873       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
[667]874         !$OMP DO
875         DO  i = nxl, nxr
876            DO  j = nys, nyn
[709]877               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
878!
879!--               Vertical heat flux
[667]880                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
881                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
882                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
883                           * w(k,j,i) * rmask(j,i,sr)
884                  IF ( humidity )  THEN
885                     pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
886                                q(k+1,j,i) - hom(k+1,1,41,sr) )
887                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
888                                                      rmask(j,i,sr)
889                  ENDIF
890               ENDDO
891            ENDDO
892         ENDDO
893
894       ENDIF
895
[1]896!
[97]897!--    Density at top follows Neumann condition
[388]898       IF ( ocean )  THEN
899          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
900          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
901       ENDIF
[97]902
903!
[1]904!--    Divergence of vertical flux of resolved scale energy and pressure
[106]905!--    fluctuations as well as flux of pressure fluctuation itself (68).
906!--    First calculate the products, then the divergence.
[1]907!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
[106]908       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
[1]909
910          sums_ll = 0.0  ! local array
911
912          !$OMP DO
913          DO  i = nxl, nxr
914             DO  j = nys, nyn
[132]915                DO  k = nzb_s_inner(j,i)+1, nzt
[1]916
917                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
918                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
[678]919                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
[1]920                           ) )**2                                          &
921                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
[678]922                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
[1]923                           ) )**2                                          &
924                   + w(k,j,i)**2                                  )
925
926                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
927                                               * ( p(k,j,i) + p(k+1,j,i) )
928
929                ENDDO
930             ENDDO
931          ENDDO
932          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
933          sums_ll(nzt+1,1) = 0.0
934          sums_ll(0,2)     = 0.0
935          sums_ll(nzt+1,2) = 0.0
936
[678]937          DO  k = nzb+1, nzt
[1]938             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
939             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
[106]940             sums_l(k,68,tn) = sums_ll(k,2)
[1]941          ENDDO
942          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
943          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
[106]944          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
[1]945
946       ENDIF
947
948!
[106]949!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
950       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
[1]951
952          !$OMP DO
953          DO  i = nxl, nxr
954             DO  j = nys, nyn
[132]955                DO  k = nzb_s_inner(j,i)+1, nzt
[1]956
[106]957                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
[1]958                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
959                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
[106]960                                                             ) * ddzw(k)
[1]961
[106]962                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
963                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
964                                                              )
965
[1]966                ENDDO
967             ENDDO
968          ENDDO
969          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
[106]970          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
[1]971
972       ENDIF
973
974!
975!--    Horizontal heat fluxes (subgrid, resolved, total).
976!--    Do it only, if profiles shall be plotted.
977       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
978
979          !$OMP DO
980          DO  i = nxl, nxr
981             DO  j = nys, nyn
[132]982                DO  k = nzb_s_inner(j,i)+1, nzt
[1]983!
984!--                Subgrid horizontal heat fluxes u"pt", v"pt"
985                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
986                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
987                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
988                                                 * ddx * rmask(j,i,sr)
989                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
990                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
991                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
992                                                 * ddy * rmask(j,i,sr)
993!
994!--                Resolved horizontal heat fluxes u*pt*, v*pt*
995                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
996                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
997                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
998                                                 pt(k,j,i)   - hom(k,1,4,sr) )
999                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1000                                 pt(k,j,i)   - hom(k,1,4,sr) )
1001                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
1002                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
1003                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1004                                                 pt(k,j,i)   - hom(k,1,4,sr) )
1005                ENDDO
1006             ENDDO
1007          ENDDO
1008!
1009!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
[97]1010          sums_l(nzb,58,tn) = 0.0
1011          sums_l(nzb,59,tn) = 0.0
1012          sums_l(nzb,60,tn) = 0.0
1013          sums_l(nzb,61,tn) = 0.0
1014          sums_l(nzb,62,tn) = 0.0
1015          sums_l(nzb,63,tn) = 0.0
[1]1016
1017       ENDIF
[87]1018
1019!
1020!--    Calculate the user-defined profiles
1021       CALL user_statistics( 'profiles', sr, tn )
[1]1022       !$OMP END PARALLEL
1023
1024!
1025!--    Summation of thread sums
1026       IF ( threads_per_task > 1 )  THEN
1027          DO  i = 1, threads_per_task-1
1028             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1029             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
[87]1030             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1031                                      sums_l(:,45:pr_palm,i)
1032             IF ( max_pr_user > 0 )  THEN
1033                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1034                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1035                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1036             ENDIF
[1]1037          ENDDO
1038       ENDIF
1039
1040#if defined( __parallel )
[667]1041
[1]1042!
1043!--    Compute total sum from local sums
[622]1044       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]1045       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
1046                           MPI_SUM, comm2d, ierr )
1047#else
1048       sums = sums_l(:,:,0)
1049#endif
1050
1051!
1052!--    Final values are obtained by division by the total number of grid points
1053!--    used for summation. After that store profiles.
1054!--    Profiles:
1055       DO  k = nzb, nzt+1
[132]1056          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
[142]1057          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
[132]1058          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
1059          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
1060          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
[142]1061          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
1062          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
[132]1063          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
1064          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
1065          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
1066          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
1067          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
1068          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
1069          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
[1]1070       ENDDO
[667]1071
[1]1072!--    Upstream-parts
[87]1073       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
[1]1074!--    u* and so on
[87]1075!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
[1]1076!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1077!--    above the topography, they are being divided by ngp_2dh(sr)
[87]1078       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
[1]1079                                    ngp_2dh(sr)
[197]1080       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1081                                    ngp_2dh(sr)
[1]1082!--    eges, e*
[87]1083       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
[132]1084                                    ngp_3d(sr)
[1]1085!--    Old and new divergence
[87]1086       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
[1]1087                                    ngp_3d_inner(sr)
1088
[87]1089!--    User-defined profiles
1090       IF ( max_pr_user > 0 )  THEN
1091          DO  k = nzb, nzt+1
1092             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1093                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
[132]1094                                    ngp_2dh_s_inner(k,sr)
[87]1095          ENDDO
1096       ENDIF
[1007]1097
[1]1098!
1099!--    Collect horizontal average in hom.
1100!--    Compute deduced averages (e.g. total heat flux)
1101       hom(:,1,3,sr)  = sums(:,3)      ! w
1102       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1103       hom(:,1,9,sr)  = sums(:,9)      ! km
1104       hom(:,1,10,sr) = sums(:,10)     ! kh
1105       hom(:,1,11,sr) = sums(:,11)     ! l
1106       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1107       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1108       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1109       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1110       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1111       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1112       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1113       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1114       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1115       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1116       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
[96]1117                                       ! profile 24 is initial profile (sa)
1118                                       ! profiles 25-29 left empty for initial
[1]1119                                       ! profiles
1120       hom(:,1,30,sr) = sums(:,30)     ! u*2
1121       hom(:,1,31,sr) = sums(:,31)     ! v*2
1122       hom(:,1,32,sr) = sums(:,32)     ! w*2
1123       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1124       hom(:,1,34,sr) = sums(:,34)     ! e*
1125       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1126       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1127       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1128       hom(:,1,38,sr) = sums(:,38)     ! w*3
[699]1129       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
[1]1130       hom(:,1,40,sr) = sums(:,40)     ! p
[531]1131       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
[1]1132       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1133       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1134       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1135       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1136       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1137       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1138       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1139       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1140       hom(:,1,54,sr) = sums(:,54)     ! ql
1141       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1142       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
[106]1143       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
[1]1144       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1145       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1146       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1147       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1148       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1149       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
[96]1150       hom(:,1,64,sr) = sums(:,64)     ! rho
1151       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1152       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1153       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
[106]1154       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1155       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
[197]1156       hom(:,1,70,sr) = sums(:,70)     ! q*2
[388]1157       hom(:,1,71,sr) = sums(:,71)     ! prho
[531]1158       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
[1053]1159       hom(:,1,73,sr) = sums(:,73)     ! nr
1160       hom(:,1,74,sr) = sums(:,74)     ! qr
1161       hom(:,1,75,sr) = sums(:,75)     ! qc
1162       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
[1179]1163                                       ! 77 is initial density profile
[1241]1164       hom(:,1,78,sr) = ug             ! ug
1165       hom(:,1,79,sr) = vg             ! vg
[1]1166
[87]1167       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
[1]1168                                       ! upstream-parts u_x, u_y, u_z, v_x,
1169                                       ! v_y, usw. (in last but one profile)
[667]1170       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
[1]1171                                       ! u*, w'u', w'v', t* (in last profile)
1172
[87]1173       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1174          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1175                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1176       ENDIF
1177
[1]1178!
1179!--    Determine the boundary layer height using two different schemes.
[94]1180!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1181!--    first relative minimum (maximum) of the total heat flux.
1182!--    The corresponding height is assumed as the boundary layer height, if it
1183!--    is less than 1.5 times the height where the heat flux becomes negative
1184!--    (positive) for the first time.
[1]1185       z_i(1) = 0.0
1186       first = .TRUE.
[667]1187
[97]1188       IF ( ocean )  THEN
1189          DO  k = nzt, nzb+1, -1
[667]1190             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1191                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
[97]1192                first = .FALSE.
1193                height = zw(k)
1194             ENDIF
1195             IF ( hom(k,1,18,sr) < 0.0  .AND. &
[667]1196                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
[97]1197                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1198                IF ( zw(k) < 1.5 * height )  THEN
1199                   z_i(1) = zw(k)
1200                ELSE
1201                   z_i(1) = height
1202                ENDIF
1203                EXIT
1204             ENDIF
1205          ENDDO
1206       ELSE
[94]1207          DO  k = nzb, nzt-1
[667]1208             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1209                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
[94]1210                first = .FALSE.
1211                height = zw(k)
[1]1212             ENDIF
[667]1213             IF ( hom(k,1,18,sr) < 0.0  .AND. & 
1214                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
[94]1215                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1216                IF ( zw(k) < 1.5 * height )  THEN
1217                   z_i(1) = zw(k)
1218                ELSE
1219                   z_i(1) = height
1220                ENDIF
1221                EXIT
1222             ENDIF
1223          ENDDO
[97]1224       ENDIF
[1]1225
1226!
[291]1227!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1228!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1229!--    maximal local temperature gradient: starting from the second (the last
1230!--    but one) vertical gridpoint, the local gradient must be at least
1231!--    0.2K/100m and greater than the next four gradients.
1232!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1233!--             ocean case!
[1]1234       z_i(2) = 0.0
[291]1235       DO  k = nzb+1, nzt+1
1236          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1237       ENDDO
1238       dptdz_threshold = 0.2 / 100.0
1239
[97]1240       IF ( ocean )  THEN
[291]1241          DO  k = nzt+1, nzb+5, -1
1242             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1243                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1244                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1245                z_i(2) = zw(k-1)
[97]1246                EXIT
1247             ENDIF
1248          ENDDO
1249       ELSE
[291]1250          DO  k = nzb+1, nzt-3
1251             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1252                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1253                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1254                z_i(2) = zw(k-1)
[97]1255                EXIT
1256             ENDIF
1257          ENDDO
1258       ENDIF
[1]1259
[87]1260       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1261       hom(nzb+7,1,pr_palm,sr) = z_i(2)
[1]1262
1263!
1264!--    Computation of both the characteristic vertical velocity and
1265!--    the characteristic convective boundary layer temperature.
1266!--    The horizontal average at nzb+1 is input for the average temperature.
[667]1267       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
1268           .AND.  z_i(1) /= 0.0 )  THEN
[87]1269          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
[94]1270                                       hom(nzb,1,18,sr) *      &
1271                                       ABS( z_i(1) ) )**0.333333333
[1]1272!--       so far this only works if Prandtl layer is used
[87]1273          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
[1]1274       ELSE
[87]1275          hom(nzb+8,1,pr_palm,sr)  = 0.0
1276          hom(nzb+11,1,pr_palm,sr) = 0.0
[1]1277       ENDIF
1278
[48]1279!
1280!--    Collect the time series quantities
[87]1281       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1282       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
[48]1283       ts_value(3,sr) = dt_3d
[87]1284       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1285       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
[48]1286       ts_value(6,sr) = u_max
1287       ts_value(7,sr) = v_max
1288       ts_value(8,sr) = w_max
[87]1289       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1290       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1291       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1292       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1293       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
[48]1294       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1295       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1296       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1297       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1298       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
[197]1299       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1300       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
[343]1301       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
[197]1302
[48]1303       IF ( ts_value(5,sr) /= 0.0 )  THEN
1304          ts_value(22,sr) = ts_value(4,sr)**2 / &
1305                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1306       ELSE
1307          ts_value(22,sr) = 10000.0
1308       ENDIF
[1]1309
[343]1310       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
[1]1311!
[48]1312!--    Calculate additional statistics provided by the user interface
[87]1313       CALL user_statistics( 'time_series', sr, 0 )
[1]1314
[48]1315    ENDDO    ! loop of the subregions
1316
[1]1317!
1318!-- If required, sum up horizontal averages for subsequent time averaging
1319    IF ( do_sum )  THEN
1320       IF ( average_count_pr == 0 )  hom_sum = 0.0
1321       hom_sum = hom_sum + hom(:,1,:,:)
1322       average_count_pr = average_count_pr + 1
1323       do_sum = .FALSE.
1324    ENDIF
1325
1326!
1327!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1328!-- This flag is reset after each time step in time_integration.
1329    flow_statistics_called = .TRUE.
1330
1331    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1332
1333
1334 END SUBROUTINE flow_statistics
[1221]1335
1336
1337#else
1338
1339
1340!------------------------------------------------------------------------------!
1341! flow statistics - accelerator version
1342!------------------------------------------------------------------------------!
1343 SUBROUTINE flow_statistics
1344
1345    USE arrays_3d
1346    USE cloud_parameters
1347    USE control_parameters
1348    USE cpulog
1349    USE grid_variables
1350    USE indices
1351    USE interfaces
1352    USE pegrid
1353    USE statistics
1354
1355    IMPLICIT NONE
1356
1357    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
1358    LOGICAL ::  first
1359    REAL    ::  dptdz_threshold, height, pts, sums_l_eper, sums_l_etot, ust, &
1360                ust2, u2, vst, vst2, v2, w2, z_i(2)
1361    REAL    ::  s1, s2, s3, s4, s5, s6, s7
1362    REAL    ::  dptdz(nzb+1:nzt+1)
1363    REAL    ::  sums_ll(nzb:nzt+1,2)
1364
1365    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
1366
1367!
1368!-- To be on the safe side, check whether flow_statistics has already been
1369!-- called once after the current time step
1370    IF ( flow_statistics_called )  THEN
1371
1372       message_string = 'flow_statistics is called two times within one ' // &
1373                        'timestep'
1374       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
1375
1376    ENDIF
1377
1378    !$acc data copyin( hom ) create( sums, sums_l )
1379
1380!
1381!-- Compute statistics for each (sub-)region
1382    DO  sr = 0, statistic_regions
1383
1384!
1385!--    Initialize (local) summation array
1386       sums_l = 0.0
1387
1388!
1389!--    Store sums that have been computed in other subroutines in summation
1390!--    array
1391       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
1392!--    WARNING: next line still has to be adjusted for OpenMP
1393       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
1394       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
1395       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
1396
1397!
1398!--    Copy the turbulent quantities, evaluated in the advection routines to
1399!--    the local array sums_l() for further computations
1400       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
1401
1402!
1403!--       According to the Neumann bc for the horizontal velocity components,
1404!--       the corresponding fluxes has to satisfiy the same bc.
1405          IF ( ocean )  THEN
1406             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
1407             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
1408          ENDIF
1409
1410          DO  i = 0, threads_per_task-1
1411!
1412!--          Swap the turbulent quantities evaluated in advec_ws.
1413             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
1414             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
1415             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
1416             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
1417             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
1418             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        &
1419                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
1420                                sums_ws2_ws_l(:,i) )    ! e*
1421             DO  k = nzb, nzt
1422                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
1423                                                      sums_us2_ws_l(k,i) +  &
1424                                                      sums_vs2_ws_l(k,i) +  &
1425                                                      sums_ws2_ws_l(k,i) )
1426             ENDDO
1427          ENDDO
1428
1429       ENDIF
1430
1431       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
1432
1433          DO  i = 0, threads_per_task-1
1434             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
1435             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
1436             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
1437                                                   sums_wsqs_ws_l(:,i) !w*q*
1438          ENDDO
1439
1440       ENDIF
1441!
1442!--    Horizontally averaged profiles of horizontal velocities and temperature.
1443!--    They must have been computed before, because they are already required
1444!--    for other horizontal averages.
1445       tn = 0
1446
1447       !$OMP PARALLEL PRIVATE( i, j, k, tn )
1448#if defined( __intel_openmp_bug )
1449       tn = omp_get_thread_num()
1450#else
1451!$     tn = omp_get_thread_num()
1452#endif
1453
1454       !$acc update device( sums_l )
1455
1456       !$OMP DO
1457       !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
1458       DO  k = nzb, nzt+1
1459          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1460          DO  i = nxl, nxr
1461             DO  j =  nys, nyn
1462!
1463!--             k+1 is used in rflags since rflags is set 0 at surface points
1464                s1 = s1 + u(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1465                s2 = s2 + v(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1466                s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1467             ENDDO
1468          ENDDO
1469          sums_l(k,1,tn) = s1
1470          sums_l(k,2,tn) = s2
1471          sums_l(k,4,tn) = s3
1472       ENDDO
[1257]1473       !$acc end parallel loop
[1221]1474
1475!
1476!--    Horizontally averaged profile of salinity
1477       IF ( ocean )  THEN
1478          !$OMP DO
1479          !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
1480          DO  k = nzb, nzt+1
1481             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1482             DO  i = nxl, nxr
1483                DO  j =  nys, nyn
1484                   s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1485                ENDDO
1486             ENDDO
1487             sums_l(k,23,tn) = s1
1488          ENDDO
[1257]1489          !$acc end parallel loop
[1221]1490       ENDIF
1491
1492!
1493!--    Horizontally averaged profiles of virtual potential temperature,
1494!--    total water content, specific humidity and liquid water potential
1495!--    temperature
1496       IF ( humidity )  THEN
1497
1498          !$OMP DO
1499          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
1500          DO  k = nzb, nzt+1
1501             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1502             DO  i = nxl, nxr
1503                DO  j =  nys, nyn
1504                   s1 = s1 + q(k,j,i)   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1505                   s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1506                ENDDO
1507             ENDDO
1508             sums_l(k,41,tn) = s1
1509             sums_l(k,44,tn) = s2
1510          ENDDO
[1257]1511          !$acc end parallel loop
[1221]1512
1513          IF ( cloud_physics )  THEN
1514             !$OMP DO
1515             !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
1516             DO  k = nzb, nzt+1
1517                !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1518                DO  i = nxl, nxr
1519                   DO  j =  nys, nyn
1520                      s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
1521                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1522                      s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * &
1523                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1524                   ENDDO
1525                ENDDO
1526                sums_l(k,42,tn) = s1
1527                sums_l(k,43,tn) = s2
1528             ENDDO
[1257]1529             !$acc end parallel loop
[1221]1530          ENDIF
1531       ENDIF
1532
1533!
1534!--    Horizontally averaged profiles of passive scalar
1535       IF ( passive_scalar )  THEN
1536          !$OMP DO
1537          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
1538          DO  k = nzb, nzt+1
1539             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1540             DO  i = nxl, nxr
1541                DO  j =  nys, nyn
1542                   s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1543                ENDDO
1544             ENDDO
1545             sums_l(k,41,tn) = s1
1546          ENDDO
[1257]1547          !$acc end parallel loop
[1221]1548       ENDIF
1549       !$OMP END PARALLEL
1550
1551!
1552!--    Summation of thread sums
1553       IF ( threads_per_task > 1 )  THEN
1554          DO  i = 1, threads_per_task-1
1555             !$acc parallel present( sums_l )
1556             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
1557             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
1558             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
1559             !$acc end parallel
1560             IF ( ocean )  THEN
1561                !$acc parallel present( sums_l )
1562                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
1563                !$acc end parallel
1564             ENDIF
1565             IF ( humidity )  THEN
1566                !$acc parallel present( sums_l )
1567                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1568                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
1569                !$acc end parallel
1570                IF ( cloud_physics )  THEN
1571                   !$acc parallel present( sums_l )
1572                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
1573                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
1574                   !$acc end parallel
1575                ENDIF
1576             ENDIF
1577             IF ( passive_scalar )  THEN
1578                !$acc parallel present( sums_l )
1579                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1580                !$acc end parallel
1581             ENDIF
1582          ENDDO
1583       ENDIF
1584
1585#if defined( __parallel )
1586!
1587!--    Compute total sum from local sums
1588       !$acc update host( sums_l )
1589       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1590       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
1591                           MPI_SUM, comm2d, ierr )
1592       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1593       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
1594                           MPI_SUM, comm2d, ierr )
1595       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1596       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
1597                           MPI_SUM, comm2d, ierr )
1598       IF ( ocean )  THEN
1599          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1600          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
1601                              MPI_REAL, MPI_SUM, comm2d, ierr )
1602       ENDIF
1603       IF ( humidity ) THEN
1604          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1605          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
1606                              MPI_REAL, MPI_SUM, comm2d, ierr )
1607          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1608          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
1609                              MPI_REAL, MPI_SUM, comm2d, ierr )
1610          IF ( cloud_physics ) THEN
1611             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1612             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
1613                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1614             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1615             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
1616                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1617          ENDIF
1618       ENDIF
1619
1620       IF ( passive_scalar )  THEN
1621          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1622          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
1623                              MPI_REAL, MPI_SUM, comm2d, ierr )
1624       ENDIF
1625       !$acc update device( sums )
1626#else
1627       !$acc parallel present( sums, sums_l )
1628       sums(:,1) = sums_l(:,1,0)
1629       sums(:,2) = sums_l(:,2,0)
1630       sums(:,4) = sums_l(:,4,0)
1631       !$acc end parallel
1632       IF ( ocean )  THEN
1633          !$acc parallel present( sums, sums_l )
1634          sums(:,23) = sums_l(:,23,0)
1635          !$acc end parallel
1636       ENDIF
1637       IF ( humidity )  THEN
1638          !$acc parallel present( sums, sums_l )
1639          sums(:,44) = sums_l(:,44,0)
1640          sums(:,41) = sums_l(:,41,0)
1641          !$acc end parallel
1642          IF ( cloud_physics )  THEN
1643             !$acc parallel present( sums, sums_l )
1644             sums(:,42) = sums_l(:,42,0)
1645             sums(:,43) = sums_l(:,43,0)
1646             !$acc end parallel
1647          ENDIF
1648       ENDIF
1649       IF ( passive_scalar )  THEN
1650          !$acc parallel present( sums, sums_l )
1651          sums(:,41) = sums_l(:,41,0)
1652          !$acc end parallel
1653       ENDIF
1654#endif
1655
1656!
1657!--    Final values are obtained by division by the total number of grid points
1658!--    used for summation. After that store profiles.
1659       !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
1660       sums(:,1) = sums(:,1) / ngp_2dh(sr)
1661       sums(:,2) = sums(:,2) / ngp_2dh(sr)
1662       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
1663       hom(:,1,1,sr) = sums(:,1)             ! u
1664       hom(:,1,2,sr) = sums(:,2)             ! v
1665       hom(:,1,4,sr) = sums(:,4)             ! pt
1666       !$acc end parallel
1667
1668!
1669!--    Salinity
1670       IF ( ocean )  THEN
1671          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1672          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
1673          hom(:,1,23,sr) = sums(:,23)             ! sa
1674          !$acc end parallel
1675       ENDIF
1676
1677!
1678!--    Humidity and cloud parameters
1679       IF ( humidity ) THEN
1680          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1681          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
1682          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1683          hom(:,1,44,sr) = sums(:,44)                ! vpt
1684          hom(:,1,41,sr) = sums(:,41)                ! qv (q)
1685          !$acc end parallel
1686          IF ( cloud_physics ) THEN
1687             !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1688             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
1689             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
1690             hom(:,1,42,sr) = sums(:,42)             ! qv
1691             hom(:,1,43,sr) = sums(:,43)             ! pt
1692             !$acc end parallel
1693          ENDIF
1694       ENDIF
1695
1696!
1697!--    Passive scalar
1698       IF ( passive_scalar )  THEN
1699          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1700          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1701          hom(:,1,41,sr) = sums(:,41)                ! s (q)
1702          !$acc end parallel
1703       ENDIF
1704
1705!
1706!--    Horizontally averaged profiles of the remaining prognostic variables,
1707!--    variances, the total and the perturbation energy (single values in last
1708!--    column of sums_l) and some diagnostic quantities.
1709!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
1710!--    ----  speaking the following k-loop would have to be split up and
1711!--          rearranged according to the staggered grid.
1712!--          However, this implies no error since staggered velocity components
1713!--          are zero at the walls and inside buildings.
1714       tn = 0
1715#if defined( __intel_openmp_bug )
1716       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
1717       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
1718       tn = omp_get_thread_num()
1719#else
1720       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
1721!$     tn = omp_get_thread_num()
1722#endif
1723       !$OMP DO
1724       !$acc parallel loop gang present( e, hom, kh, km, p, pt, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, s5, s6, s7 )
1725       DO  k = nzb, nzt+1
1726          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
1727          DO  i = nxl, nxr
1728             DO  j =  nys, nyn
1729!
1730!--             Prognostic and diagnostic variables
1731                s1 = s1 + w(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1732                s2 = s2 + e(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1733                s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1734                s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1735                s5 = s5 + p(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1736                s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * &
1737                          rflags_invers(j,i,k+1)
1738!
1739!--             Higher moments
1740!--             (Computation of the skewness of w further below)
1741                s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1742             ENDDO
1743          ENDDO
1744          sums_l(k,3,tn)  = s1
1745          sums_l(k,8,tn)  = s2
1746          sums_l(k,9,tn)  = s3
1747          sums_l(k,10,tn) = s4
1748          sums_l(k,40,tn) = s5
1749          sums_l(k,33,tn) = s6
1750          sums_l(k,38,tn) = s7
1751       ENDDO
[1257]1752       !$acc end parallel loop
[1221]1753
1754       IF ( humidity )  THEN
1755          !$OMP DO
1756          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
1757          DO  k = nzb, nzt+1
1758             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1759             DO  i = nxl, nxr
1760                DO  j =  nys, nyn
1761                   s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &
1762                             rflags_invers(j,i,k+1)
1763                ENDDO
1764             ENDDO
1765             sums_l(k,70,tn) = s1
1766          ENDDO
[1257]1767          !$acc end parallel loop
[1221]1768       ENDIF
1769
1770!
1771!--    Total and perturbation energy for the total domain (being
1772!--    collected in the last column of sums_l).
1773       !$OMP DO
1774       !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
1775       DO  i = nxl, nxr
1776          DO  j =  nys, nyn
1777             DO  k = nzb, nzt+1
1778                s1 = s1 + 0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) * &
1779                          rmask(j,i,sr) * rflags_invers(j,i,k+1)
1780             ENDDO
1781          ENDDO
1782       ENDDO
[1257]1783       !$acc end parallel loop
[1221]1784       !$acc parallel present( sums_l )
1785       sums_l(nzb+4,pr_palm,tn) = s1
1786       !$acc end parallel
1787
1788       !$OMP DO
1789       !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
1790       !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
1791       DO  i = nxl, nxr
1792          DO  j =  nys, nyn
1793!
1794!--          2D-arrays (being collected in the last column of sums_l)
1795             s1 = s1 + us(j,i)   * rmask(j,i,sr)
1796             s2 = s2 + usws(j,i) * rmask(j,i,sr)
1797             s3 = s3 + vsws(j,i) * rmask(j,i,sr)
1798             s4 = s4 + ts(j,i)   * rmask(j,i,sr)
1799          ENDDO
1800       ENDDO
1801       sums_l(nzb,pr_palm,tn)   = s1
1802       sums_l(nzb+1,pr_palm,tn) = s2
1803       sums_l(nzb+2,pr_palm,tn) = s3
1804       sums_l(nzb+3,pr_palm,tn) = s4
1805       !$acc end parallel
1806
1807       IF ( humidity )  THEN
1808          !$acc parallel present( qs, rmask, sums_l ) create( s1 )
1809          !$acc loop vector collapse( 2 ) reduction( +: s1 )
1810          DO  i = nxl, nxr
1811             DO  j =  nys, nyn
1812                s1 = s1 + qs(j,i) * rmask(j,i,sr)
1813             ENDDO
1814          ENDDO
1815          sums_l(nzb+12,pr_palm,tn) = s1
1816          !$acc end parallel
1817       ENDIF
1818
1819!
1820!--    Computation of statistics when ws-scheme is not used. Else these
1821!--    quantities are evaluated in the advection routines.
1822       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0 )  THEN
1823
1824          !$OMP DO
1825          !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
1826          DO  k = nzb, nzt+1
1827             !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
1828             DO  i = nxl, nxr
1829                DO  j =  nys, nyn
1830                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
1831                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
1832                   w2   = w(k,j,i)**2
1833
1834                   s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1835                   s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1836                   s3 = s3 + w2   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1837!
1838!--                Perturbation energy
1839                   s4 = s4 + 0.5 * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * &
1840                             rflags_invers(j,i,k+1)
1841                ENDDO
1842             ENDDO
1843             sums_l(k,30,tn) = s1
1844             sums_l(k,31,tn) = s2
1845             sums_l(k,32,tn) = s3
1846             sums_l(k,34,tn) = s4
1847          ENDDO
[1257]1848          !$acc end parallel loop
[1221]1849!
1850!--       Total perturbation TKE
1851          !$OMP DO
1852          !$acc parallel present( sums_l ) create( s1 )
1853          !$acc loop reduction( +: s1 )
1854          DO  k = nzb, nzt+1
1855             s1 = s1 + sums_l(k,34,tn)
1856          ENDDO
1857          sums_l(nzb+5,pr_palm,tn) = s1
1858          !$acc end parallel
1859
1860       ENDIF
1861
1862!
1863!--    Horizontally averaged profiles of the vertical fluxes
1864
1865!
1866!--    Subgridscale fluxes.
1867!--    WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
1868!--    -------  should be calculated there in a different way. This is done
1869!--             in the next loop further below, where results from this loop are
1870!--             overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
1871!--             The non-flat case still has to be handled.
1872!--    NOTE: for simplicity, nzb_s_inner is used below, although
1873!--    ----  strictly speaking the following k-loop would have to be
1874!--          split up according to the staggered grid.
1875!--          However, this implies no error since staggered velocity
1876!--          components are zero at the walls and inside buildings.
1877       !$OMP DO
1878       !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
1879       DO  k = nzb, nzt_diff
1880          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1881          DO  i = nxl, nxr
1882             DO  j = nys, nyn
1883
1884!
1885!--             Momentum flux w"u"
1886                s1 = s1 - 0.25 * (                   &
1887                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
1888                                                           ) * (               &
1889                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
1890                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
1891                                                               )               &
1892                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1893!
1894!--             Momentum flux w"v"
1895                s2 = s2 - 0.25 * (                   &
1896                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
1897                                                           ) * (               &
1898                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
1899                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
1900                                                               )               &
1901                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1902!
1903!--             Heat flux w"pt"
1904                s3 = s3 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1905                              * ( pt(k+1,j,i) - pt(k,j,i) )   &
1906                              * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1907             ENDDO
1908          ENDDO
1909          sums_l(k,12,tn) = s1
1910          sums_l(k,14,tn) = s2
1911          sums_l(k,16,tn) = s3
1912       ENDDO
[1257]1913       !$acc end parallel loop
[1221]1914
1915!
1916!--    Salinity flux w"sa"
1917       IF ( ocean )  THEN
1918          !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
1919          DO  k = nzb, nzt_diff
1920             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1921             DO  i = nxl, nxr
1922                DO  j = nys, nyn
1923                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1924                                 * ( sa(k+1,j,i) - sa(k,j,i) )   &
1925                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1926                ENDDO
1927             ENDDO
1928             sums_l(k,65,tn) = s1
1929          ENDDO
[1257]1930          !$acc end parallel loop
[1221]1931       ENDIF
1932
1933!
1934!--    Buoyancy flux, water flux (humidity flux) w"q"
1935       IF ( humidity ) THEN
1936
1937          !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
1938          DO  k = nzb, nzt_diff
1939             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1940             DO  i = nxl, nxr
1941                DO  j = nys, nyn
1942                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1943                                 * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
1944                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1945                   s2 = s2 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1946                                 * ( q(k+1,j,i) - q(k,j,i) )     &
1947                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1948                ENDDO
1949             ENDDO
1950             sums_l(k,45,tn) = s1
1951             sums_l(k,48,tn) = s2
1952          ENDDO
[1257]1953          !$acc end parallel loop
[1221]1954
1955          IF ( cloud_physics ) THEN
1956
1957             !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
1958             DO  k = nzb, nzt_diff
1959                !$acc loop vector collapse( 2 ) reduction( +: s1 )
1960                DO  i = nxl, nxr
1961                   DO  j = nys, nyn
1962                      s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )    &
1963                                    * ( ( q(k+1,j,i) - ql(k+1,j,i) ) &
1964                                      - ( q(k,j,i) - ql(k,j,i) ) )   &
1965                                    * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1966                   ENDDO
1967                ENDDO
1968                sums_l(k,51,tn) = s1
1969             ENDDO
[1257]1970             !$acc end parallel loop
[1221]1971
1972          ENDIF
1973
1974       ENDIF
1975!
1976!--    Passive scalar flux
1977       IF ( passive_scalar )  THEN
1978
1979          !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
1980          DO  k = nzb, nzt_diff
1981             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1982             DO  i = nxl, nxr
1983                DO  j = nys, nyn
1984                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1985                                 * ( q(k+1,j,i) - q(k,j,i) )     &
1986                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1987                ENDDO
1988             ENDDO
1989             sums_l(k,48,tn) = s1
1990          ENDDO
[1257]1991          !$acc end parallel loop
[1221]1992
1993       ENDIF
1994
1995       IF ( use_surface_fluxes )  THEN
1996
1997          !$OMP DO
1998          !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
1999          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2000          DO  i = nxl, nxr
2001             DO  j =  nys, nyn
2002!
2003!--             Subgridscale fluxes in the Prandtl layer
2004                s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
2005                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
2006                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
2007                s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
2008                s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
2009             ENDDO
2010          ENDDO
2011          sums_l(nzb,12,tn) = s1
2012          sums_l(nzb,14,tn) = s2
2013          sums_l(nzb,16,tn) = s3
2014          sums_l(nzb,58,tn) = s4
2015          sums_l(nzb,61,tn) = s5
2016          !$acc end parallel
2017
2018          IF ( ocean )  THEN
2019
2020             !$OMP DO
2021             !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
2022             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2023             DO  i = nxl, nxr
2024                DO  j =  nys, nyn
2025                   s1 = s1 + saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
2026                ENDDO
2027             ENDDO
2028             sums_l(nzb,65,tn) = s1
2029             !$acc end parallel
2030
2031          ENDIF
2032
2033          IF ( humidity )  THEN
2034
2035             !$OMP DO
2036             !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
2037             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2038             DO  i = nxl, nxr
2039                DO  j =  nys, nyn
2040                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2041                   s2 = s2 + ( ( 1.0 + 0.61 * q(nzb,j,i) ) * shf(j,i) &
2042                               + 0.61 * pt(nzb,j,i) * qsws(j,i) )
2043                ENDDO
2044             ENDDO
2045             sums_l(nzb,48,tn) = s1
2046             sums_l(nzb,45,tn) = s2
2047             !$acc end parallel
2048
2049             IF ( cloud_droplets )  THEN
2050
2051                !$OMP DO
2052                !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
2053                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2054                DO  i = nxl, nxr
2055                   DO  j =  nys, nyn
2056                      s1 = s1 + ( ( 1.0 + 0.61 * q(nzb,j,i) - ql(nzb,j,i) ) * &
2057                                  shf(j,i) + 0.61 * pt(nzb,j,i) * qsws(j,i) )
2058                   ENDDO
2059                ENDDO
2060                sums_l(nzb,45,tn) = s1
2061                !$acc end parallel
2062
2063             ENDIF
2064
2065             IF ( cloud_physics )  THEN
2066
2067                !$OMP DO
2068                !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2069                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2070                DO  i = nxl, nxr
2071                   DO  j =  nys, nyn
2072!
2073!--                   Formula does not work if ql(nzb) /= 0.0
2074                      s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
2075                   ENDDO
2076                ENDDO
2077                sums_l(nzb,51,tn) = s1
2078                !$acc end parallel
2079
2080             ENDIF
2081
2082          ENDIF
2083
2084          IF ( passive_scalar )  THEN
2085
2086             !$OMP DO
2087             !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2088             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2089             DO  i = nxl, nxr
2090                DO  j =  nys, nyn
2091                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2092                ENDDO
2093             ENDDO
2094             sums_l(nzb,48,tn) = s1
2095             !$acc end parallel
2096
2097          ENDIF
2098
2099       ENDIF
2100
2101!
2102!--    Subgridscale fluxes at the top surface
2103       IF ( use_top_fluxes )  THEN
2104
2105          !$OMP DO
2106          !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
2107          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2108          DO  i = nxl, nxr
2109             DO  j =  nys, nyn
2110                s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
2111                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
2112                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
2113                s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
2114                s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
2115             ENDDO
2116          ENDDO
2117          sums_l(nzt:nzt+1,12,tn) = s1
2118          sums_l(nzt:nzt+1,14,tn) = s2
2119          sums_l(nzt:nzt+1,16,tn) = s3
2120          sums_l(nzt:nzt+1,58,tn) = s4
2121          sums_l(nzt:nzt+1,61,tn) = s5
2122          !$acc end parallel
2123
2124          IF ( ocean )  THEN
2125
2126             !$OMP DO
2127             !$acc parallel present( rmask, saswst, sums_l ) create( s1 )
2128             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2129             DO  i = nxl, nxr
2130                DO  j =  nys, nyn
2131                   s1 = s1 + saswst(j,i) * rmask(j,i,sr)  ! w"sa"
2132                ENDDO
2133             ENDDO
2134             sums_l(nzt,65,tn) = s1
2135             !$acc end parallel
2136
2137          ENDIF
2138
2139          IF ( humidity )  THEN
2140
2141             !$OMP DO
2142             !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
2143             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2144             DO  i = nxl, nxr
2145                DO  j =  nys, nyn
2146                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2147                   s2 = s2 + ( ( 1.0 + 0.61 * q(nzt,j,i) ) * tswst(j,i) + &
2148                               0.61 * pt(nzt,j,i) * qswst(j,i) )
2149                ENDDO
2150             ENDDO
2151             sums_l(nzt,48,tn) = s1
2152             sums_l(nzt,45,tn) = s2
2153             !$acc end parallel
2154
2155             IF ( cloud_droplets )  THEN
2156
2157                !$OMP DO
2158                !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
2159                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2160                DO  i = nxl, nxr
2161                   DO  j =  nys, nyn
2162                      s1 = s1 + ( ( 1.0 + 0.61 * q(nzt,j,i) - ql(nzt,j,i) ) * &
2163                                  tswst(j,i) + 0.61 * pt(nzt,j,i) * qswst(j,i) )
2164                   ENDDO
2165                ENDDO
2166                sums_l(nzt,45,tn) = s1
2167                !$acc end parallel
2168
2169             ENDIF
2170
2171             IF ( cloud_physics )  THEN
2172
2173                !$OMP DO
2174                !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2175                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2176                DO  i = nxl, nxr
2177                   DO  j =  nys, nyn
2178!
2179!--                   Formula does not work if ql(nzb) /= 0.0
2180                      s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2181                   ENDDO
2182                ENDDO
2183                sums_l(nzt,51,tn) = s1
2184                !$acc end parallel
2185
2186             ENDIF
2187
2188          ENDIF
2189
2190          IF ( passive_scalar )  THEN
2191
2192             !$OMP DO
2193             !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2194             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2195             DO  i = nxl, nxr
2196                DO  j =  nys, nyn
2197                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2198                ENDDO
2199             ENDDO
2200             sums_l(nzt,48,tn) = s1
2201             !$acc end parallel
2202
2203          ENDIF
2204
2205       ENDIF
2206
2207!
2208!--    Resolved fluxes (can be computed for all horizontal points)
2209!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
2210!--    ----  speaking the following k-loop would have to be split up and
2211!--          rearranged according to the staggered grid.
2212       !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
2213       DO  k = nzb, nzt_diff
2214          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2215          DO  i = nxl, nxr
2216             DO  j = nys, nyn
2217                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
2218                              u(k+1,j,i) - hom(k+1,1,1,sr) )
2219                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
2220                              v(k+1,j,i) - hom(k+1,1,2,sr) )
2221                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
2222                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
2223!
2224!--             Higher moments
2225                s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2226                s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2227!
2228!--             Energy flux w*e* (has to be adjusted?)
2229                s3 = s3 + w(k,j,i) * 0.5 * ( ust**2 + vst**2 + w(k,j,i)**2 ) &
2230                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2231             ENDDO
2232          ENDDO
2233          sums_l(k,35,tn) = s1
2234          sums_l(k,36,tn) = s2
2235          sums_l(k,37,tn) = s3
2236       ENDDO
[1257]2237       !$acc end parallel loop
[1221]2238
2239!
2240!--    Salinity flux and density (density does not belong to here,
2241!--    but so far there is no other suitable place to calculate)
2242       IF ( ocean )  THEN
2243
2244          IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
2245
2246             !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
2247             DO  k = nzb, nzt_diff
2248                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2249                DO  i = nxl, nxr
2250                   DO  j = nys, nyn
2251                      s1 = s1 + 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) +   &
2252                                        sa(k+1,j,i) - hom(k+1,1,23,sr) ) &
2253                                    * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2254                   ENDDO
2255                ENDDO
2256                sums_l(k,66,tn) = s1
2257             ENDDO
[1257]2258             !$acc end parallel loop
[1221]2259
2260          ENDIF
2261
2262          !$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
2263          DO  k = nzb, nzt_diff
2264             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2265             DO  i = nxl, nxr
2266                DO  j = nys, nyn
2267                   s1 = s1 + rho(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2268                   s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2269                ENDDO
2270             ENDDO
2271             sums_l(k,64,tn) = s1
2272             sums_l(k,71,tn) = s2
2273          ENDDO
[1257]2274          !$acc end parallel loop
[1221]2275
2276       ENDIF
2277
2278!
2279!--    Buoyancy flux, water flux, humidity flux, liquid water
2280!--    content, rain drop concentration and rain water content
2281       IF ( humidity )  THEN
2282
2283          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
2284
2285             !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2286             DO  k = nzb, nzt_diff
2287                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2288                DO  i = nxl, nxr
2289                   DO  j = nys, nyn
2290                      s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
2291                                        vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
2292                                      w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2293                   ENDDO
2294                ENDDO
2295                sums_l(k,46,tn) = s1
2296             ENDDO
[1257]2297             !$acc end parallel loop
[1221]2298
2299             IF ( .NOT. cloud_droplets )  THEN
2300
2301                !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
2302                DO  k = nzb, nzt_diff
2303                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2304                   DO  i = nxl, nxr
2305                      DO  j = nys, nyn
2306                         s1 = s1 + 0.5 * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
2307                                           ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
2308                                       * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2309                      ENDDO
2310                   ENDDO
2311                   sums_l(k,52,tn) = s1
2312                ENDDO
[1257]2313                !$acc end parallel loop
[1221]2314
2315                IF ( icloud_scheme == 0  )  THEN
2316
2317                   !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
2318                   DO  k = nzb, nzt_diff
2319                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2320                      DO  i = nxl, nxr
2321                         DO  j = nys, nyn
2322                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2323                            s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2324                         ENDDO
2325                      ENDDO
2326                      sums_l(k,54,tn) = s1
2327                      sums_l(k,75,tn) = s2
2328                   ENDDO
[1257]2329                   !$acc end parallel loop
[1221]2330
2331                   IF ( precipitation )  THEN
2332
2333                      !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2334                      DO  k = nzb, nzt_diff
2335                         !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2336                         DO  i = nxl, nxr
2337                            DO  j = nys, nyn
2338                               s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2339                               s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2340                               s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2341                            ENDDO
2342                         ENDDO
2343                         sums_l(k,73,tn) = s1
2344                         sums_l(k,74,tn) = s2
2345                         sums_l(k,76,tn) = s3
2346                      ENDDO
[1257]2347                      !$acc end parallel loop
[1221]2348
2349                   ENDIF
2350
2351                ELSE
2352
2353                   !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2354                   DO  k = nzb, nzt_diff
2355                      !$acc loop vector collapse( 2 ) reduction( +: s1 )
2356                      DO  i = nxl, nxr
2357                         DO  j = nys, nyn
2358                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2359                         ENDDO
2360                      ENDDO
2361                      sums_l(k,54,tn) = s1
2362                   ENDDO
[1257]2363                   !$acc end parallel loop
[1221]2364
2365                ENDIF
2366
2367             ELSE
2368
2369                !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2370                DO  k = nzb, nzt_diff
2371                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2372                   DO  i = nxl, nxr
2373                      DO  j = nys, nyn
2374                         s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2375                      ENDDO
2376                   ENDDO
2377                   sums_l(k,54,tn) = s1
2378                ENDDO
[1257]2379                !$acc end parallel loop
[1221]2380
2381             ENDIF
2382
2383          ELSE
2384
2385             IF( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2386
2387                !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2388                DO  k = nzb, nzt_diff
2389                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2390                   DO  i = nxl, nxr
2391                      DO  j = nys, nyn
2392                         s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
2393                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
2394                                       * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2395                      ENDDO
2396                   ENDDO
2397                   sums_l(k,46,tn) = s1
2398                ENDDO
[1257]2399                !$acc end parallel loop
[1221]2400
2401             ELSEIF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
2402
2403                !$acc parallel loop present( hom, sums_l )
2404                DO  k = nzb, nzt_diff
2405                   sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
2406                                             0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
2407                ENDDO
[1257]2408                !$acc end parallel loop
[1221]2409
2410             ENDIF
2411
2412          ENDIF
2413
2414       ENDIF
2415!
2416!--    Passive scalar flux
2417       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
2418
2419          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2420          DO  k = nzb, nzt_diff
2421             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2422             DO  i = nxl, nxr
2423                DO  j = nys, nyn
2424                   s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
2425                                     q(k+1,j,i) - hom(k+1,1,41,sr) ) &
2426                                 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2427                ENDDO
2428             ENDDO
2429             sums_l(k,49,tn) = s1
2430          ENDDO
[1257]2431          !$acc end parallel loop
[1221]2432
2433       ENDIF
2434
2435!
2436!--    For speed optimization fluxes which have been computed in part directly
2437!--    inside the WS advection routines are treated seperatly
2438!--    Momentum fluxes first:
2439       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
2440
2441          !$OMP DO
2442          !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2443          DO  k = nzb, nzt_diff
2444             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2445             DO  i = nxl, nxr
2446                DO  j = nys, nyn
2447                   ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
2448                               u(k+1,j,i) - hom(k+1,1,1,sr) )
2449                   vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
2450                               v(k+1,j,i) - hom(k+1,1,2,sr) )
2451!
2452!--                Momentum flux w*u*
2453                   s1 = s1 + 0.5 * ( w(k,j,i-1) + w(k,j,i) ) &
2454                                 * ust * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2455!
2456!--                Momentum flux w*v*
2457                   s2 = s2 + 0.5 * ( w(k,j-1,i) + w(k,j,i) ) &
2458                                 * vst * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2459                ENDDO
2460             ENDDO
2461             sums_l(k,13,tn) = s1
2462             sums_l(k,15,tn) = s1
2463          ENDDO
[1257]2464          !$acc end parallel loop
[1221]2465
2466       ENDIF
2467
2468       IF ( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2469
2470          !$OMP DO
2471          !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
2472          DO  k = nzb, nzt_diff
2473             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2474             DO  i = nxl, nxr
2475                DO  j = nys, nyn
2476!
2477!--                Vertical heat flux
2478                   s1 = s1 + 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
2479                                     pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
2480                                 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2481                ENDDO
2482             ENDDO
2483             sums_l(k,17,tn) = s1
2484          ENDDO
[1257]2485          !$acc end parallel loop
[1221]2486
2487          IF ( humidity )  THEN
2488
2489             !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2490             DO  k = nzb, nzt_diff
2491                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2492                DO  i = nxl, nxr
2493                   DO  j = nys, nyn
2494                      s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
2495                                        q(k+1,j,i) - hom(k+1,1,41,sr) ) &
2496                                    * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2497                   ENDDO
2498                ENDDO
2499                sums_l(k,49,tn) = s1
2500             ENDDO
[1257]2501             !$acc end parallel loop
[1221]2502
2503          ENDIF
2504
2505       ENDIF
2506
2507
2508!
2509!--    Density at top follows Neumann condition
2510       IF ( ocean )  THEN
2511          !$acc parallel present( sums_l )
2512          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
2513          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
2514          !$acc end parallel
2515       ENDIF
2516
2517!
2518!--    Divergence of vertical flux of resolved scale energy and pressure
2519!--    fluctuations as well as flux of pressure fluctuation itself (68).
2520!--    First calculate the products, then the divergence.
2521!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
2522       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
2523
2524          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
2525          sums_ll = 0.0  ! local array
2526
2527          !$OMP DO
2528          DO  i = nxl, nxr
2529             DO  j = nys, nyn
2530                DO  k = nzb_s_inner(j,i)+1, nzt
2531
2532                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
2533                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
2534                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
2535                           ) )**2                                          &
2536                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
2537                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
2538                           ) )**2                                          &
2539                   + w(k,j,i)**2                                  )
2540
2541                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
2542                                               * ( p(k,j,i) + p(k+1,j,i) )
2543
2544                ENDDO
2545             ENDDO
2546          ENDDO
2547          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
2548          sums_ll(nzt+1,1) = 0.0
2549          sums_ll(0,2)     = 0.0
2550          sums_ll(nzt+1,2) = 0.0
2551
2552          DO  k = nzb+1, nzt
2553             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
2554             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
2555             sums_l(k,68,tn) = sums_ll(k,2)
2556          ENDDO
2557          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
2558          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
2559          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
2560
2561       ENDIF
2562
2563!
2564!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
2565       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
2566
2567          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
2568          !$OMP DO
2569          DO  i = nxl, nxr
2570             DO  j = nys, nyn
2571                DO  k = nzb_s_inner(j,i)+1, nzt
2572
2573                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
2574                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2575                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
2576                                                             ) * ddzw(k)
2577
2578                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
2579                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2580                                                              )
2581
2582                ENDDO
2583             ENDDO
2584          ENDDO
2585          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
2586          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
2587
2588       ENDIF
2589
2590!
2591!--    Horizontal heat fluxes (subgrid, resolved, total).
2592!--    Do it only, if profiles shall be plotted.
2593       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
2594
2595          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
2596          !$OMP DO
2597          DO  i = nxl, nxr
2598             DO  j = nys, nyn
2599                DO  k = nzb_s_inner(j,i)+1, nzt
2600!
2601!--                Subgrid horizontal heat fluxes u"pt", v"pt"
2602                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
2603                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
2604                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
2605                                                 * ddx * rmask(j,i,sr)
2606                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
2607                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
2608                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
2609                                                 * ddy * rmask(j,i,sr)
2610!
2611!--                Resolved horizontal heat fluxes u*pt*, v*pt*
2612                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
2613                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
2614                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
2615                                                 pt(k,j,i)   - hom(k,1,4,sr) )
2616                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
2617                                 pt(k,j,i)   - hom(k,1,4,sr) )
2618                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
2619                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
2620                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
2621                                                 pt(k,j,i)   - hom(k,1,4,sr) )
2622                ENDDO
2623             ENDDO
2624          ENDDO
2625!
2626!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
2627          sums_l(nzb,58,tn) = 0.0
2628          sums_l(nzb,59,tn) = 0.0
2629          sums_l(nzb,60,tn) = 0.0
2630          sums_l(nzb,61,tn) = 0.0
2631          sums_l(nzb,62,tn) = 0.0
2632          sums_l(nzb,63,tn) = 0.0
2633
2634       ENDIF
2635
2636!
2637!--    Calculate the user-defined profiles
2638       CALL user_statistics( 'profiles', sr, tn )
2639       !$OMP END PARALLEL
2640
2641!
2642!--    Summation of thread sums
2643       IF ( threads_per_task > 1 )  THEN
2644          STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
2645          DO  i = 1, threads_per_task-1
2646             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
2647             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
2648             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
2649                                      sums_l(:,45:pr_palm,i)
2650             IF ( max_pr_user > 0 )  THEN
2651                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
2652                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
2653                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
2654             ENDIF
2655          ENDDO
2656       ENDIF
2657
2658       !$acc update host( hom, sums, sums_l )
2659
2660#if defined( __parallel )
2661
2662!
2663!--    Compute total sum from local sums
2664       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2665       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
2666                           MPI_SUM, comm2d, ierr )
2667#else
2668       sums = sums_l(:,:,0)
2669#endif
2670
2671!
2672!--    Final values are obtained by division by the total number of grid points
2673!--    used for summation. After that store profiles.
2674!--    Profiles:
2675       DO  k = nzb, nzt+1
2676          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
2677          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
2678          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
2679          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
2680          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
2681          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
2682          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
2683          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
2684          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
2685          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
2686          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
2687          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
2688          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
2689          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
2690       ENDDO
2691
2692!--    Upstream-parts
2693       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
2694!--    u* and so on
2695!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
2696!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
2697!--    above the topography, they are being divided by ngp_2dh(sr)
2698       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
2699                                    ngp_2dh(sr)
2700       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
2701                                    ngp_2dh(sr)
2702!--    eges, e*
2703       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
2704                                    ngp_3d(sr)
2705!--    Old and new divergence
2706       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
2707                                    ngp_3d_inner(sr)
2708
2709!--    User-defined profiles
2710       IF ( max_pr_user > 0 )  THEN
2711          DO  k = nzb, nzt+1
2712             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
2713                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
2714                                    ngp_2dh_s_inner(k,sr)
2715          ENDDO
2716       ENDIF
2717
2718!
2719!--    Collect horizontal average in hom.
2720!--    Compute deduced averages (e.g. total heat flux)
2721       hom(:,1,3,sr)  = sums(:,3)      ! w
2722       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
2723       hom(:,1,9,sr)  = sums(:,9)      ! km
2724       hom(:,1,10,sr) = sums(:,10)     ! kh
2725       hom(:,1,11,sr) = sums(:,11)     ! l
2726       hom(:,1,12,sr) = sums(:,12)     ! w"u"
2727       hom(:,1,13,sr) = sums(:,13)     ! w*u*
2728       hom(:,1,14,sr) = sums(:,14)     ! w"v"
2729       hom(:,1,15,sr) = sums(:,15)     ! w*v*
2730       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
2731       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
2732       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
2733       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
2734       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
2735       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
2736       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
2737                                       ! profile 24 is initial profile (sa)
2738                                       ! profiles 25-29 left empty for initial
2739                                       ! profiles
2740       hom(:,1,30,sr) = sums(:,30)     ! u*2
2741       hom(:,1,31,sr) = sums(:,31)     ! v*2
2742       hom(:,1,32,sr) = sums(:,32)     ! w*2
2743       hom(:,1,33,sr) = sums(:,33)     ! pt*2
2744       hom(:,1,34,sr) = sums(:,34)     ! e*
2745       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
2746       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
2747       hom(:,1,37,sr) = sums(:,37)     ! w*e*
2748       hom(:,1,38,sr) = sums(:,38)     ! w*3
2749       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
2750       hom(:,1,40,sr) = sums(:,40)     ! p
2751       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
2752       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
2753       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
2754       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
2755       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
2756       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
2757       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
2758       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
2759       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
2760       hom(:,1,54,sr) = sums(:,54)     ! ql
2761       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
2762       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
2763       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
2764       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
2765       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
2766       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
2767       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
2768       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
2769       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
2770       hom(:,1,64,sr) = sums(:,64)     ! rho
2771       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
2772       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
2773       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
2774       hom(:,1,68,sr) = sums(:,68)     ! w*p*
2775       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
2776       hom(:,1,70,sr) = sums(:,70)     ! q*2
2777       hom(:,1,71,sr) = sums(:,71)     ! prho
2778       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
2779       hom(:,1,73,sr) = sums(:,73)     ! nr
2780       hom(:,1,74,sr) = sums(:,74)     ! qr
2781       hom(:,1,75,sr) = sums(:,75)     ! qc
2782       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
2783                                       ! 77 is initial density profile
[1241]2784       hom(:,1,78,sr) = ug             ! ug
2785       hom(:,1,79,sr) = vg             ! vg
[1221]2786
2787       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
2788                                       ! upstream-parts u_x, u_y, u_z, v_x,
2789                                       ! v_y, usw. (in last but one profile)
2790       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
2791                                       ! u*, w'u', w'v', t* (in last profile)
2792
2793       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
2794          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
2795                               sums(:,pr_palm+1:pr_palm+max_pr_user)
2796       ENDIF
2797
2798!
2799!--    Determine the boundary layer height using two different schemes.
2800!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
2801!--    first relative minimum (maximum) of the total heat flux.
2802!--    The corresponding height is assumed as the boundary layer height, if it
2803!--    is less than 1.5 times the height where the heat flux becomes negative
2804!--    (positive) for the first time.
2805       z_i(1) = 0.0
2806       first = .TRUE.
2807
2808       IF ( ocean )  THEN
2809          DO  k = nzt, nzb+1, -1
2810             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
2811                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
2812                first = .FALSE.
2813                height = zw(k)
2814             ENDIF
2815             IF ( hom(k,1,18,sr) < 0.0  .AND. &
2816                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
2817                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
2818                IF ( zw(k) < 1.5 * height )  THEN
2819                   z_i(1) = zw(k)
2820                ELSE
2821                   z_i(1) = height
2822                ENDIF
2823                EXIT
2824             ENDIF
2825          ENDDO
2826       ELSE
2827          DO  k = nzb, nzt-1
2828             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
2829                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
2830                first = .FALSE.
2831                height = zw(k)
2832             ENDIF
2833             IF ( hom(k,1,18,sr) < 0.0  .AND. &
2834                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
2835                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
2836                IF ( zw(k) < 1.5 * height )  THEN
2837                   z_i(1) = zw(k)
2838                ELSE
2839                   z_i(1) = height
2840                ENDIF
2841                EXIT
2842             ENDIF
2843          ENDDO
2844       ENDIF
2845
2846!
2847!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
2848!--    by Uhlenbrock(2006). The boundary layer height is the height with the
2849!--    maximal local temperature gradient: starting from the second (the last
2850!--    but one) vertical gridpoint, the local gradient must be at least
2851!--    0.2K/100m and greater than the next four gradients.
2852!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
2853!--             ocean case!
2854       z_i(2) = 0.0
2855       DO  k = nzb+1, nzt+1
2856          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
2857       ENDDO
2858       dptdz_threshold = 0.2 / 100.0
2859
2860       IF ( ocean )  THEN
2861          DO  k = nzt+1, nzb+5, -1
2862             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
2863                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
2864                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
2865                z_i(2) = zw(k-1)
2866                EXIT
2867             ENDIF
2868          ENDDO
2869       ELSE
2870          DO  k = nzb+1, nzt-3
2871             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
2872                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
2873                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
2874                z_i(2) = zw(k-1)
2875                EXIT
2876             ENDIF
2877          ENDDO
2878       ENDIF
2879
2880       hom(nzb+6,1,pr_palm,sr) = z_i(1)
2881       hom(nzb+7,1,pr_palm,sr) = z_i(2)
2882
2883!
2884!--    Computation of both the characteristic vertical velocity and
2885!--    the characteristic convective boundary layer temperature.
2886!--    The horizontal average at nzb+1 is input for the average temperature.
2887       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
2888           .AND.  z_i(1) /= 0.0 )  THEN
2889          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
2890                                       hom(nzb,1,18,sr) *      &
2891                                       ABS( z_i(1) ) )**0.333333333
2892!--       so far this only works if Prandtl layer is used
2893          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
2894       ELSE
2895          hom(nzb+8,1,pr_palm,sr)  = 0.0
2896          hom(nzb+11,1,pr_palm,sr) = 0.0
2897       ENDIF
2898
2899!
2900!--    Collect the time series quantities
2901       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
2902       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
2903       ts_value(3,sr) = dt_3d
2904       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
2905       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
2906       ts_value(6,sr) = u_max
2907       ts_value(7,sr) = v_max
2908       ts_value(8,sr) = w_max
2909       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
2910       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
2911       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
2912       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
2913       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
2914       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
2915       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
2916       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
2917       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
2918       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
2919       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
2920       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
2921       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
2922
2923       IF ( ts_value(5,sr) /= 0.0 )  THEN
2924          ts_value(22,sr) = ts_value(4,sr)**2 / &
2925                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
2926       ELSE
2927          ts_value(22,sr) = 10000.0
2928       ENDIF
2929
2930       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
2931!
2932!--    Calculate additional statistics provided by the user interface
2933       CALL user_statistics( 'time_series', sr, 0 )
2934
2935    ENDDO    ! loop of the subregions
2936
2937    !$acc end data
2938
2939!
2940!-- If required, sum up horizontal averages for subsequent time averaging
2941    IF ( do_sum )  THEN
2942       IF ( average_count_pr == 0 )  hom_sum = 0.0
2943       hom_sum = hom_sum + hom(:,1,:,:)
2944       average_count_pr = average_count_pr + 1
2945       do_sum = .FALSE.
2946    ENDIF
2947
2948!
2949!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
2950!-- This flag is reset after each time step in time_integration.
2951    flow_statistics_called = .TRUE.
2952
2953    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
2954
2955
2956 END SUBROUTINE flow_statistics
2957#endif
Note: See TracBrowser for help on using the repository browser.