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

Last change on this file since 1299 was 1299, checked in by heinze, 7 years ago

enable usage of large_scale subsidence in combination with large_scale_forcing

  • Property svn:keywords set to Id
File size: 120.7 KB
Line 
1#if ! defined( __openacc )
2 SUBROUTINE flow_statistics
3
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!
21! Current revisions:
22! -----------------
23! Output of large scale vertical velocity w_subs
24!
25! Former revisions:
26! -----------------
27! $Id: flow_statistics.f90 1299 2014-03-06 13:15:21Z heinze $
28!
29! 1257 2013-11-08 15:18:40Z raasch
30! openacc "end parallel" replaced by "end parallel loop"
31!
32! 1241 2013-10-30 11:36:58Z heinze
33! Output of ug and vg
34!
35! 1221 2013-09-10 08:59:13Z raasch
36! ported for openACC in separate #else branch
37!
38! 1179 2013-06-14 05:57:58Z raasch
39! comment for profile 77 added
40!
41! 1115 2013-03-26 18:16:16Z hoffmann
42! ql is calculated by calc_liquid_water_content
43!
44! 1111 2013-03-08 23:54:10Z raasch
45! openACC directive added
46!
47! 1053 2012-11-13 17:11:03Z hoffmann
48! additions for two-moment cloud physics scheme:
49! +nr, qr, qc, prr
50!
51! 1036 2012-10-22 13:43:42Z raasch
52! code put under GPL (PALM 3.9)
53!
54! 1007 2012-09-19 14:30:36Z franke
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
59!
60! 801 2012-01-10 17:30:36Z suehring
61! Calculation of turbulent fluxes in advec_ws is now thread-safe.
62!
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!
67! 709 2011-03-30 09:31:40Z raasch
68! formatting adjustments
69!
70! 699 2011-03-22 17:52:22Z suehring
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.
76!, tend
77! 696 2011-03-18 07:03:49Z raasch
78! Bugfix: Summation of Wicker-Skamarock scheme fluxes and variances for all
79! threads
80!
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!
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!
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!
100! 624 2010-12-10 11:46:30Z heinze
101! Calculation of q*2 added
102!
103! 622 2010-12-10 08:08:13Z raasch
104! optional barriers included in order to speed up collective operations
105!
106! 388 2009-09-23 09:40:33Z raasch
107! Vertical profiles of potential density and hydrostatic pressure are
108! calculated.
109! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
110! end.
111! Temperature gradient criterion for estimating the boundary layer height
112! replaced by the gradient criterion of Sullivan et al. (1998).
113! Output of messages replaced by message handling routine.
114!
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!
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!
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!
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!
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!
141! 82 2007-04-16 15:40:52Z raasch
142! Cpp-directive lcmuk changed to intel_openmp_bug
143!
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!
149! 19 2007-02-23 04:53:48Z raasch
150! fluxes at top modified (tswst, qswst)
151!
152! RCS Log replace by Id keyword, revision history cleaned up
153!
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!
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.
171!------------------------------------------------------------------------------!
172
173    USE arrays_3d
174    USE cloud_parameters
175    USE control_parameters
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
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)
190    REAL    ::  sums_ll(nzb:nzt+1,2)
191
192    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
193
194    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, usws, v, vsws, w )
195
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
200
201       message_string = 'flow_statistics is called two times within one ' // &
202                        'timestep'
203       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
204
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
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
223
224!
225!--    Copy the turbulent quantities, evaluated in the advection routines to
226!--    the local array sums_l() for further computations
227       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
228
229!
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
233             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
234             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
235          ENDIF
236
237          DO  i = 0, threads_per_task-1
238!
239!--          Swap the turbulent quantities evaluated in advec_ws.
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
245             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        & 
246                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
247                                sums_ws2_ws_l(:,i) )    ! e*
248             DO  k = nzb, nzt
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) )
253             ENDDO
254          ENDDO
255
256       ENDIF
257
258       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
259
260          DO  i = 0, threads_per_task-1
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*
263             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
264                                                   sums_wsqs_ws_l(:,i) !w*q*
265          ENDDO
266
267       ENDIF
268!
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
273
274       !$OMP PARALLEL PRIVATE( i, j, k, tn )
275#if defined( __intel_openmp_bug )
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
284             DO  k = nzb_s_inner(j,i), nzt+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!
293!--    Horizontally averaged profile of salinity
294       IF ( ocean )  THEN
295          !$OMP DO
296          DO  i = nxl, nxr
297             DO  j =  nys, nyn
298                DO  k = nzb_s_inner(j,i), nzt+1
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!
307!--    Horizontally averaged profiles of virtual potential temperature,
308!--    total water content, specific humidity and liquid water potential
309!--    temperature
310       IF ( humidity )  THEN
311          !$OMP DO
312          DO  i = nxl, nxr
313             DO  j =  nys, nyn
314                DO  k = nzb_s_inner(j,i), nzt+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
326                   DO  k = nzb_s_inner(j,i), nzt+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
344                DO  k = nzb_s_inner(j,i), nzt+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)
358             IF ( ocean )  THEN
359                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
360             ENDIF
361             IF ( humidity )  THEN
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
378       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
379       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
380                           MPI_SUM, comm2d, ierr )
381       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
382       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
383                           MPI_SUM, comm2d, ierr )
384       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
385       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
386                           MPI_SUM, comm2d, ierr )
387       IF ( ocean )  THEN
388          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
389          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
390                              MPI_REAL, MPI_SUM, comm2d, ierr )
391       ENDIF
392       IF ( humidity ) THEN
393          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
394          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
395                              MPI_REAL, MPI_SUM, comm2d, ierr )
396          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
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
400             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
401             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
402                                 MPI_REAL, MPI_SUM, comm2d, ierr )
403             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
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
410          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
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)
418       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
419       IF ( humidity ) THEN
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.
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)
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
440
441!
442!--    Salinity
443       IF ( ocean )  THEN
444          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
445          hom(:,1,23,sr) = sums(:,23)             ! sa
446       ENDIF
447
448!
449!--    Humidity and cloud parameters
450       IF ( humidity ) THEN
451          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
452          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
453          hom(:,1,44,sr) = sums(:,44)             ! vpt
454          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
455          IF ( cloud_physics ) THEN
456             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
457             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
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
465       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /  &
466            ngp_2dh_s_inner(:,sr)                    ! s (q)
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.
472!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
473!--    ----  speaking the following k-loop would have to be split up and
474!--          rearranged according to the staggered grid.
475!--          However, this implies no error since staggered velocity components
476!--          are zero at the walls and inside buildings.
477       tn = 0
478#if defined( __intel_openmp_bug )
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
490             DO  k = nzb_s_inner(j,i), nzt+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)
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
506
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)
511
512                sums_l_etot  = sums_l_etot + &
513                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
514                                        w(k,j,i)**2 ) * rmask(j,i,sr)
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.
521             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
522!
523!--          2D-arrays (being collected in the last column of sums_l)
524             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +   &
525                                        us(j,i)   * rmask(j,i,sr)
526             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
527                                        usws(j,i) * rmask(j,i,sr)
528             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
529                                        vsws(j,i) * rmask(j,i,sr)
530             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
531                                        ts(j,i)   * rmask(j,i,sr)
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
536          ENDDO
537       ENDDO
538
539!
540!--    Computation of statistics when ws-scheme is not used. Else these
541!--    quantities are evaluated in the advection routines.
542       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 )  THEN
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
571
572!
573!--    Horizontally averaged profiles of the vertical fluxes
574
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)
581!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
582!--          ----  strictly speaking the following k-loop would have to be
583!--                split up according to the staggered grid.
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
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!
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!
622!--             Buoyancy flux, water flux (humidity flux) w"q"
623                IF ( humidity ) THEN
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)
632
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"
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
670                IF ( humidity )  THEN
671                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
672                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
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
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!
697!--          Subgridscale fluxes at the top surface
698             IF ( use_top_fluxes )  THEN
699                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
700                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
701                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
702                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
703                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
704                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
705                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
706                                    0.0 * rmask(j,i,sr)           ! u"pt"
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
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
714                IF ( humidity )  THEN
715                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
716                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
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
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) + &
736                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
737                ENDIF
738             ENDIF
739
740!
741!--          Resolved fluxes (can be computed for all horizontal points)
742!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
743!--          ----  speaking the following k-loop would have to be split up and
744!--                rearranged according to the staggered grid.
745             DO  k = nzb_s_inner(j,i), nzt
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) )
752
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!
760!--             Salinity flux and density (density does not belong to here,
761!--             but so far there is no other suitable place to calculate)
762                IF ( ocean )  THEN
763                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
764                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
765                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
766                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
767                                                       rmask(j,i,sr)
768                   ENDIF
769                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
770                                                       rmask(j,i,sr)
771                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
772                                                       rmask(j,i,sr)
773                ENDIF
774
775!
776!--             Buoyancy flux, water flux, humidity flux, liquid water
777!--             content, rain drop concentration and rain water content
778                IF ( humidity )  THEN
779                   IF ( cloud_physics .OR. cloud_droplets )  THEN
780                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +        &
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) * &
783                                                          rmask(j,i,sr)
784                      IF ( .NOT. cloud_droplets )  THEN
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) ) -                   &
789                              hom(k+1,1,42,sr) )
790                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
791                                                             rmask(j,i,sr)
792                         IF ( icloud_scheme == 0  )  THEN
793                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
794                                                                rmask(j,i,sr)
795                            sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *    &
796                                                                rmask(j,i,sr)
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
805                         ELSE
806                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
807                                                                rmask(j,i,sr)
808                         ENDIF
809                      ELSE
810                         sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *       &
811                                                             rmask(j,i,sr)
812                      ENDIF
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
820                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) *  &
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
825                ENDIF
826!
827!--             Passive scalar flux
828                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca   &
829                     .OR. sr /= 0 ) )  THEN
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*
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)
842             ENDDO
843          ENDDO
844       ENDDO
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:
849       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
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) )
858!
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
871
872       ENDIF
873       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
874         !$OMP DO
875         DO  i = nxl, nxr
876            DO  j = nys, nyn
877               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
878!
879!--               Vertical heat flux
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
896!
897!--    Density at top follows Neumann condition
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
902
903!
904!--    Divergence of vertical flux of resolved scale energy and pressure
905!--    fluctuations as well as flux of pressure fluctuation itself (68).
906!--    First calculate the products, then the divergence.
907!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
908       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
909
910          sums_ll = 0.0  ! local array
911
912          !$OMP DO
913          DO  i = nxl, nxr
914             DO  j = nys, nyn
915                DO  k = nzb_s_inner(j,i)+1, nzt
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)   &
919                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
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)   &
922                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
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
937          DO  k = nzb+1, nzt
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)
940             sums_l(k,68,tn) = sums_ll(k,2)
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)
944          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
945
946       ENDIF
947
948!
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
951
952          !$OMP DO
953          DO  i = nxl, nxr
954             DO  j = nys, nyn
955                DO  k = nzb_s_inner(j,i)+1, nzt
956
957                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
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)   &
960                                                             ) * ddzw(k)
961
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
966                ENDDO
967             ENDDO
968          ENDDO
969          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
970          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
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
982                DO  k = nzb_s_inner(j,i)+1, nzt
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)
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
1016
1017       ENDIF
1018
1019!
1020!--    Calculate the user-defined profiles
1021       CALL user_statistics( 'profiles', sr, tn )
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)
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
1037          ENDDO
1038       ENDIF
1039
1040#if defined( __parallel )
1041
1042!
1043!--    Compute total sum from local sums
1044       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
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
1056          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
1057          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
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)
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)
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)
1070       ENDDO
1071
1072!--    Upstream-parts
1073       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
1074!--    u* and so on
1075!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
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)
1078       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
1079                                    ngp_2dh(sr)
1080       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1081                                    ngp_2dh(sr)
1082!--    eges, e*
1083       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
1084                                    ngp_3d(sr)
1085!--    Old and new divergence
1086       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1087                                    ngp_3d_inner(sr)
1088
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) / &
1094                                    ngp_2dh_s_inner(k,sr)
1095          ENDDO
1096       ENDIF
1097
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
1117                                       ! profile 24 is initial profile (sa)
1118                                       ! profiles 25-29 left empty for initial
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
1129       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
1130       hom(:,1,40,sr) = sums(:,40)     ! p
1131       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
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
1143       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
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
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
1154       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1155       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
1156       hom(:,1,70,sr) = sums(:,70)     ! q*2
1157       hom(:,1,71,sr) = sums(:,71)     ! prho
1158       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
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)
1163                                       ! 77 is initial density profile
1164       hom(:,1,78,sr) = ug             ! ug
1165       hom(:,1,79,sr) = vg             ! vg
1166       hom(:,1,80,sr) = w_subs         ! w_subs
1167
1168       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
1169                                       ! upstream-parts u_x, u_y, u_z, v_x,
1170                                       ! v_y, usw. (in last but one profile)
1171       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1172                                       ! u*, w'u', w'v', t* (in last profile)
1173
1174       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1175          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1176                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1177       ENDIF
1178
1179!
1180!--    Determine the boundary layer height using two different schemes.
1181!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1182!--    first relative minimum (maximum) of the total heat flux.
1183!--    The corresponding height is assumed as the boundary layer height, if it
1184!--    is less than 1.5 times the height where the heat flux becomes negative
1185!--    (positive) for the first time.
1186       z_i(1) = 0.0
1187       first = .TRUE.
1188
1189       IF ( ocean )  THEN
1190          DO  k = nzt, nzb+1, -1
1191             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1192                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
1193                first = .FALSE.
1194                height = zw(k)
1195             ENDIF
1196             IF ( hom(k,1,18,sr) < 0.0  .AND. &
1197                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1198                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1199                IF ( zw(k) < 1.5 * height )  THEN
1200                   z_i(1) = zw(k)
1201                ELSE
1202                   z_i(1) = height
1203                ENDIF
1204                EXIT
1205             ENDIF
1206          ENDDO
1207       ELSE
1208          DO  k = nzb, nzt-1
1209             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1210                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
1211                first = .FALSE.
1212                height = zw(k)
1213             ENDIF
1214             IF ( hom(k,1,18,sr) < 0.0  .AND. & 
1215                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1216                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1217                IF ( zw(k) < 1.5 * height )  THEN
1218                   z_i(1) = zw(k)
1219                ELSE
1220                   z_i(1) = height
1221                ENDIF
1222                EXIT
1223             ENDIF
1224          ENDDO
1225       ENDIF
1226
1227!
1228!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1229!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1230!--    maximal local temperature gradient: starting from the second (the last
1231!--    but one) vertical gridpoint, the local gradient must be at least
1232!--    0.2K/100m and greater than the next four gradients.
1233!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1234!--             ocean case!
1235       z_i(2) = 0.0
1236       DO  k = nzb+1, nzt+1
1237          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1238       ENDDO
1239       dptdz_threshold = 0.2 / 100.0
1240
1241       IF ( ocean )  THEN
1242          DO  k = nzt+1, nzb+5, -1
1243             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1244                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1245                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1246                z_i(2) = zw(k-1)
1247                EXIT
1248             ENDIF
1249          ENDDO
1250       ELSE
1251          DO  k = nzb+1, nzt-3
1252             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1253                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1254                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1255                z_i(2) = zw(k-1)
1256                EXIT
1257             ENDIF
1258          ENDDO
1259       ENDIF
1260
1261       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1262       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1263
1264!
1265!--    Computation of both the characteristic vertical velocity and
1266!--    the characteristic convective boundary layer temperature.
1267!--    The horizontal average at nzb+1 is input for the average temperature.
1268       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
1269           .AND.  z_i(1) /= 0.0 )  THEN
1270          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
1271                                       hom(nzb,1,18,sr) *      &
1272                                       ABS( z_i(1) ) )**0.333333333
1273!--       so far this only works if Prandtl layer is used
1274          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
1275       ELSE
1276          hom(nzb+8,1,pr_palm,sr)  = 0.0
1277          hom(nzb+11,1,pr_palm,sr) = 0.0
1278       ENDIF
1279
1280!
1281!--    Collect the time series quantities
1282       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1283       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
1284       ts_value(3,sr) = dt_3d
1285       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1286       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
1287       ts_value(6,sr) = u_max
1288       ts_value(7,sr) = v_max
1289       ts_value(8,sr) = w_max
1290       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1291       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1292       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1293       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1294       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
1295       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1296       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1297       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1298       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1299       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
1300       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1301       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
1302       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1303
1304       IF ( ts_value(5,sr) /= 0.0 )  THEN
1305          ts_value(22,sr) = ts_value(4,sr)**2 / &
1306                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1307       ELSE
1308          ts_value(22,sr) = 10000.0
1309       ENDIF
1310
1311       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1312!
1313!--    Calculate additional statistics provided by the user interface
1314       CALL user_statistics( 'time_series', sr, 0 )
1315
1316    ENDDO    ! loop of the subregions
1317
1318!
1319!-- If required, sum up horizontal averages for subsequent time averaging
1320    IF ( do_sum )  THEN
1321       IF ( average_count_pr == 0 )  hom_sum = 0.0
1322       hom_sum = hom_sum + hom(:,1,:,:)
1323       average_count_pr = average_count_pr + 1
1324       do_sum = .FALSE.
1325    ENDIF
1326
1327!
1328!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1329!-- This flag is reset after each time step in time_integration.
1330    flow_statistics_called = .TRUE.
1331
1332    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1333
1334
1335 END SUBROUTINE flow_statistics
1336
1337
1338#else
1339
1340
1341!------------------------------------------------------------------------------!
1342! flow statistics - accelerator version
1343!------------------------------------------------------------------------------!
1344 SUBROUTINE flow_statistics
1345
1346    USE arrays_3d
1347    USE cloud_parameters
1348    USE control_parameters
1349    USE cpulog
1350    USE grid_variables
1351    USE indices
1352    USE interfaces
1353    USE pegrid
1354    USE statistics
1355
1356    IMPLICIT NONE
1357
1358    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
1359    LOGICAL ::  first
1360    REAL    ::  dptdz_threshold, height, pts, sums_l_eper, sums_l_etot, ust, &
1361                ust2, u2, vst, vst2, v2, w2, z_i(2)
1362    REAL    ::  s1, s2, s3, s4, s5, s6, s7
1363    REAL    ::  dptdz(nzb+1:nzt+1)
1364    REAL    ::  sums_ll(nzb:nzt+1,2)
1365
1366    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
1367
1368!
1369!-- To be on the safe side, check whether flow_statistics has already been
1370!-- called once after the current time step
1371    IF ( flow_statistics_called )  THEN
1372
1373       message_string = 'flow_statistics is called two times within one ' // &
1374                        'timestep'
1375       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
1376
1377    ENDIF
1378
1379    !$acc data copyin( hom ) create( sums, sums_l )
1380
1381!
1382!-- Compute statistics for each (sub-)region
1383    DO  sr = 0, statistic_regions
1384
1385!
1386!--    Initialize (local) summation array
1387       sums_l = 0.0
1388
1389!
1390!--    Store sums that have been computed in other subroutines in summation
1391!--    array
1392       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
1393!--    WARNING: next line still has to be adjusted for OpenMP
1394       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
1395       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
1396       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
1397
1398!
1399!--    Copy the turbulent quantities, evaluated in the advection routines to
1400!--    the local array sums_l() for further computations
1401       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
1402
1403!
1404!--       According to the Neumann bc for the horizontal velocity components,
1405!--       the corresponding fluxes has to satisfiy the same bc.
1406          IF ( ocean )  THEN
1407             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
1408             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
1409          ENDIF
1410
1411          DO  i = 0, threads_per_task-1
1412!
1413!--          Swap the turbulent quantities evaluated in advec_ws.
1414             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
1415             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
1416             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
1417             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
1418             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
1419             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        &
1420                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
1421                                sums_ws2_ws_l(:,i) )    ! e*
1422             DO  k = nzb, nzt
1423                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
1424                                                      sums_us2_ws_l(k,i) +  &
1425                                                      sums_vs2_ws_l(k,i) +  &
1426                                                      sums_ws2_ws_l(k,i) )
1427             ENDDO
1428          ENDDO
1429
1430       ENDIF
1431
1432       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
1433
1434          DO  i = 0, threads_per_task-1
1435             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
1436             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
1437             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
1438                                                   sums_wsqs_ws_l(:,i) !w*q*
1439          ENDDO
1440
1441       ENDIF
1442!
1443!--    Horizontally averaged profiles of horizontal velocities and temperature.
1444!--    They must have been computed before, because they are already required
1445!--    for other horizontal averages.
1446       tn = 0
1447
1448       !$OMP PARALLEL PRIVATE( i, j, k, tn )
1449#if defined( __intel_openmp_bug )
1450       tn = omp_get_thread_num()
1451#else
1452!$     tn = omp_get_thread_num()
1453#endif
1454
1455       !$acc update device( sums_l )
1456
1457       !$OMP DO
1458       !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
1459       DO  k = nzb, nzt+1
1460          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1461          DO  i = nxl, nxr
1462             DO  j =  nys, nyn
1463!
1464!--             k+1 is used in rflags since rflags is set 0 at surface points
1465                s1 = s1 + u(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1466                s2 = s2 + v(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1467                s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1468             ENDDO
1469          ENDDO
1470          sums_l(k,1,tn) = s1
1471          sums_l(k,2,tn) = s2
1472          sums_l(k,4,tn) = s3
1473       ENDDO
1474       !$acc end parallel loop
1475
1476!
1477!--    Horizontally averaged profile of salinity
1478       IF ( ocean )  THEN
1479          !$OMP DO
1480          !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
1481          DO  k = nzb, nzt+1
1482             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1483             DO  i = nxl, nxr
1484                DO  j =  nys, nyn
1485                   s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1486                ENDDO
1487             ENDDO
1488             sums_l(k,23,tn) = s1
1489          ENDDO
1490          !$acc end parallel loop
1491       ENDIF
1492
1493!
1494!--    Horizontally averaged profiles of virtual potential temperature,
1495!--    total water content, specific humidity and liquid water potential
1496!--    temperature
1497       IF ( humidity )  THEN
1498
1499          !$OMP DO
1500          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
1501          DO  k = nzb, nzt+1
1502             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1503             DO  i = nxl, nxr
1504                DO  j =  nys, nyn
1505                   s1 = s1 + q(k,j,i)   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1506                   s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1507                ENDDO
1508             ENDDO
1509             sums_l(k,41,tn) = s1
1510             sums_l(k,44,tn) = s2
1511          ENDDO
1512          !$acc end parallel loop
1513
1514          IF ( cloud_physics )  THEN
1515             !$OMP DO
1516             !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
1517             DO  k = nzb, nzt+1
1518                !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1519                DO  i = nxl, nxr
1520                   DO  j =  nys, nyn
1521                      s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
1522                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1523                      s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * &
1524                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1525                   ENDDO
1526                ENDDO
1527                sums_l(k,42,tn) = s1
1528                sums_l(k,43,tn) = s2
1529             ENDDO
1530             !$acc end parallel loop
1531          ENDIF
1532       ENDIF
1533
1534!
1535!--    Horizontally averaged profiles of passive scalar
1536       IF ( passive_scalar )  THEN
1537          !$OMP DO
1538          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
1539          DO  k = nzb, nzt+1
1540             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1541             DO  i = nxl, nxr
1542                DO  j =  nys, nyn
1543                   s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1544                ENDDO
1545             ENDDO
1546             sums_l(k,41,tn) = s1
1547          ENDDO
1548          !$acc end parallel loop
1549       ENDIF
1550       !$OMP END PARALLEL
1551
1552!
1553!--    Summation of thread sums
1554       IF ( threads_per_task > 1 )  THEN
1555          DO  i = 1, threads_per_task-1
1556             !$acc parallel present( sums_l )
1557             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
1558             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
1559             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
1560             !$acc end parallel
1561             IF ( ocean )  THEN
1562                !$acc parallel present( sums_l )
1563                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
1564                !$acc end parallel
1565             ENDIF
1566             IF ( humidity )  THEN
1567                !$acc parallel present( sums_l )
1568                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1569                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
1570                !$acc end parallel
1571                IF ( cloud_physics )  THEN
1572                   !$acc parallel present( sums_l )
1573                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
1574                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
1575                   !$acc end parallel
1576                ENDIF
1577             ENDIF
1578             IF ( passive_scalar )  THEN
1579                !$acc parallel present( sums_l )
1580                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1581                !$acc end parallel
1582             ENDIF
1583          ENDDO
1584       ENDIF
1585
1586#if defined( __parallel )
1587!
1588!--    Compute total sum from local sums
1589       !$acc update host( sums_l )
1590       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1591       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
1592                           MPI_SUM, comm2d, ierr )
1593       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1594       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
1595                           MPI_SUM, comm2d, ierr )
1596       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1597       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
1598                           MPI_SUM, comm2d, ierr )
1599       IF ( ocean )  THEN
1600          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1601          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
1602                              MPI_REAL, MPI_SUM, comm2d, ierr )
1603       ENDIF
1604       IF ( humidity ) THEN
1605          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1606          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
1607                              MPI_REAL, MPI_SUM, comm2d, ierr )
1608          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1609          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
1610                              MPI_REAL, MPI_SUM, comm2d, ierr )
1611          IF ( cloud_physics ) THEN
1612             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1613             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
1614                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1615             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1616             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
1617                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1618          ENDIF
1619       ENDIF
1620
1621       IF ( passive_scalar )  THEN
1622          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1623          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
1624                              MPI_REAL, MPI_SUM, comm2d, ierr )
1625       ENDIF
1626       !$acc update device( sums )
1627#else
1628       !$acc parallel present( sums, sums_l )
1629       sums(:,1) = sums_l(:,1,0)
1630       sums(:,2) = sums_l(:,2,0)
1631       sums(:,4) = sums_l(:,4,0)
1632       !$acc end parallel
1633       IF ( ocean )  THEN
1634          !$acc parallel present( sums, sums_l )
1635          sums(:,23) = sums_l(:,23,0)
1636          !$acc end parallel
1637       ENDIF
1638       IF ( humidity )  THEN
1639          !$acc parallel present( sums, sums_l )
1640          sums(:,44) = sums_l(:,44,0)
1641          sums(:,41) = sums_l(:,41,0)
1642          !$acc end parallel
1643          IF ( cloud_physics )  THEN
1644             !$acc parallel present( sums, sums_l )
1645             sums(:,42) = sums_l(:,42,0)
1646             sums(:,43) = sums_l(:,43,0)
1647             !$acc end parallel
1648          ENDIF
1649       ENDIF
1650       IF ( passive_scalar )  THEN
1651          !$acc parallel present( sums, sums_l )
1652          sums(:,41) = sums_l(:,41,0)
1653          !$acc end parallel
1654       ENDIF
1655#endif
1656
1657!
1658!--    Final values are obtained by division by the total number of grid points
1659!--    used for summation. After that store profiles.
1660       !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
1661       sums(:,1) = sums(:,1) / ngp_2dh(sr)
1662       sums(:,2) = sums(:,2) / ngp_2dh(sr)
1663       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
1664       hom(:,1,1,sr) = sums(:,1)             ! u
1665       hom(:,1,2,sr) = sums(:,2)             ! v
1666       hom(:,1,4,sr) = sums(:,4)             ! pt
1667       !$acc end parallel
1668
1669!
1670!--    Salinity
1671       IF ( ocean )  THEN
1672          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1673          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
1674          hom(:,1,23,sr) = sums(:,23)             ! sa
1675          !$acc end parallel
1676       ENDIF
1677
1678!
1679!--    Humidity and cloud parameters
1680       IF ( humidity ) THEN
1681          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1682          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
1683          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1684          hom(:,1,44,sr) = sums(:,44)                ! vpt
1685          hom(:,1,41,sr) = sums(:,41)                ! qv (q)
1686          !$acc end parallel
1687          IF ( cloud_physics ) THEN
1688             !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1689             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
1690             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
1691             hom(:,1,42,sr) = sums(:,42)             ! qv
1692             hom(:,1,43,sr) = sums(:,43)             ! pt
1693             !$acc end parallel
1694          ENDIF
1695       ENDIF
1696
1697!
1698!--    Passive scalar
1699       IF ( passive_scalar )  THEN
1700          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1701          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1702          hom(:,1,41,sr) = sums(:,41)                ! s (q)
1703          !$acc end parallel
1704       ENDIF
1705
1706!
1707!--    Horizontally averaged profiles of the remaining prognostic variables,
1708!--    variances, the total and the perturbation energy (single values in last
1709!--    column of sums_l) and some diagnostic quantities.
1710!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
1711!--    ----  speaking the following k-loop would have to be split up and
1712!--          rearranged according to the staggered grid.
1713!--          However, this implies no error since staggered velocity components
1714!--          are zero at the walls and inside buildings.
1715       tn = 0
1716#if defined( __intel_openmp_bug )
1717       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
1718       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
1719       tn = omp_get_thread_num()
1720#else
1721       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
1722!$     tn = omp_get_thread_num()
1723#endif
1724       !$OMP DO
1725       !$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 )
1726       DO  k = nzb, nzt+1
1727          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
1728          DO  i = nxl, nxr
1729             DO  j =  nys, nyn
1730!
1731!--             Prognostic and diagnostic variables
1732                s1 = s1 + w(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1733                s2 = s2 + e(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1734                s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1735                s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1736                s5 = s5 + p(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1737                s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * &
1738                          rflags_invers(j,i,k+1)
1739!
1740!--             Higher moments
1741!--             (Computation of the skewness of w further below)
1742                s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1743             ENDDO
1744          ENDDO
1745          sums_l(k,3,tn)  = s1
1746          sums_l(k,8,tn)  = s2
1747          sums_l(k,9,tn)  = s3
1748          sums_l(k,10,tn) = s4
1749          sums_l(k,40,tn) = s5
1750          sums_l(k,33,tn) = s6
1751          sums_l(k,38,tn) = s7
1752       ENDDO
1753       !$acc end parallel loop
1754
1755       IF ( humidity )  THEN
1756          !$OMP DO
1757          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
1758          DO  k = nzb, nzt+1
1759             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1760             DO  i = nxl, nxr
1761                DO  j =  nys, nyn
1762                   s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &
1763                             rflags_invers(j,i,k+1)
1764                ENDDO
1765             ENDDO
1766             sums_l(k,70,tn) = s1
1767          ENDDO
1768          !$acc end parallel loop
1769       ENDIF
1770
1771!
1772!--    Total and perturbation energy for the total domain (being
1773!--    collected in the last column of sums_l).
1774       !$OMP DO
1775       !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
1776       DO  i = nxl, nxr
1777          DO  j =  nys, nyn
1778             DO  k = nzb, nzt+1
1779                s1 = s1 + 0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) * &
1780                          rmask(j,i,sr) * rflags_invers(j,i,k+1)
1781             ENDDO
1782          ENDDO
1783       ENDDO
1784       !$acc end parallel loop
1785       !$acc parallel present( sums_l )
1786       sums_l(nzb+4,pr_palm,tn) = s1
1787       !$acc end parallel
1788
1789       !$OMP DO
1790       !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
1791       !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
1792       DO  i = nxl, nxr
1793          DO  j =  nys, nyn
1794!
1795!--          2D-arrays (being collected in the last column of sums_l)
1796             s1 = s1 + us(j,i)   * rmask(j,i,sr)
1797             s2 = s2 + usws(j,i) * rmask(j,i,sr)
1798             s3 = s3 + vsws(j,i) * rmask(j,i,sr)
1799             s4 = s4 + ts(j,i)   * rmask(j,i,sr)
1800          ENDDO
1801       ENDDO
1802       sums_l(nzb,pr_palm,tn)   = s1
1803       sums_l(nzb+1,pr_palm,tn) = s2
1804       sums_l(nzb+2,pr_palm,tn) = s3
1805       sums_l(nzb+3,pr_palm,tn) = s4
1806       !$acc end parallel
1807
1808       IF ( humidity )  THEN
1809          !$acc parallel present( qs, rmask, sums_l ) create( s1 )
1810          !$acc loop vector collapse( 2 ) reduction( +: s1 )
1811          DO  i = nxl, nxr
1812             DO  j =  nys, nyn
1813                s1 = s1 + qs(j,i) * rmask(j,i,sr)
1814             ENDDO
1815          ENDDO
1816          sums_l(nzb+12,pr_palm,tn) = s1
1817          !$acc end parallel
1818       ENDIF
1819
1820!
1821!--    Computation of statistics when ws-scheme is not used. Else these
1822!--    quantities are evaluated in the advection routines.
1823       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0 )  THEN
1824
1825          !$OMP DO
1826          !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
1827          DO  k = nzb, nzt+1
1828             !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
1829             DO  i = nxl, nxr
1830                DO  j =  nys, nyn
1831                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
1832                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
1833                   w2   = w(k,j,i)**2
1834
1835                   s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1836                   s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1837                   s3 = s3 + w2   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1838!
1839!--                Perturbation energy
1840                   s4 = s4 + 0.5 * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * &
1841                             rflags_invers(j,i,k+1)
1842                ENDDO
1843             ENDDO
1844             sums_l(k,30,tn) = s1
1845             sums_l(k,31,tn) = s2
1846             sums_l(k,32,tn) = s3
1847             sums_l(k,34,tn) = s4
1848          ENDDO
1849          !$acc end parallel loop
1850!
1851!--       Total perturbation TKE
1852          !$OMP DO
1853          !$acc parallel present( sums_l ) create( s1 )
1854          !$acc loop reduction( +: s1 )
1855          DO  k = nzb, nzt+1
1856             s1 = s1 + sums_l(k,34,tn)
1857          ENDDO
1858          sums_l(nzb+5,pr_palm,tn) = s1
1859          !$acc end parallel
1860
1861       ENDIF
1862
1863!
1864!--    Horizontally averaged profiles of the vertical fluxes
1865
1866!
1867!--    Subgridscale fluxes.
1868!--    WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
1869!--    -------  should be calculated there in a different way. This is done
1870!--             in the next loop further below, where results from this loop are
1871!--             overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
1872!--             The non-flat case still has to be handled.
1873!--    NOTE: for simplicity, nzb_s_inner is used below, although
1874!--    ----  strictly speaking the following k-loop would have to be
1875!--          split up according to the staggered grid.
1876!--          However, this implies no error since staggered velocity
1877!--          components are zero at the walls and inside buildings.
1878       !$OMP DO
1879       !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
1880       DO  k = nzb, nzt_diff
1881          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1882          DO  i = nxl, nxr
1883             DO  j = nys, nyn
1884
1885!
1886!--             Momentum flux w"u"
1887                s1 = s1 - 0.25 * (                   &
1888                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
1889                                                           ) * (               &
1890                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
1891                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
1892                                                               )               &
1893                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1894!
1895!--             Momentum flux w"v"
1896                s2 = s2 - 0.25 * (                   &
1897                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
1898                                                           ) * (               &
1899                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
1900                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
1901                                                               )               &
1902                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1903!
1904!--             Heat flux w"pt"
1905                s3 = s3 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1906                              * ( pt(k+1,j,i) - pt(k,j,i) )   &
1907                              * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1908             ENDDO
1909          ENDDO
1910          sums_l(k,12,tn) = s1
1911          sums_l(k,14,tn) = s2
1912          sums_l(k,16,tn) = s3
1913       ENDDO
1914       !$acc end parallel loop
1915
1916!
1917!--    Salinity flux w"sa"
1918       IF ( ocean )  THEN
1919          !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
1920          DO  k = nzb, nzt_diff
1921             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1922             DO  i = nxl, nxr
1923                DO  j = nys, nyn
1924                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1925                                 * ( sa(k+1,j,i) - sa(k,j,i) )   &
1926                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1927                ENDDO
1928             ENDDO
1929             sums_l(k,65,tn) = s1
1930          ENDDO
1931          !$acc end parallel loop
1932       ENDIF
1933
1934!
1935!--    Buoyancy flux, water flux (humidity flux) w"q"
1936       IF ( humidity ) THEN
1937
1938          !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
1939          DO  k = nzb, nzt_diff
1940             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1941             DO  i = nxl, nxr
1942                DO  j = nys, nyn
1943                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1944                                 * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
1945                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1946                   s2 = s2 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1947                                 * ( q(k+1,j,i) - q(k,j,i) )     &
1948                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1949                ENDDO
1950             ENDDO
1951             sums_l(k,45,tn) = s1
1952             sums_l(k,48,tn) = s2
1953          ENDDO
1954          !$acc end parallel loop
1955
1956          IF ( cloud_physics ) THEN
1957
1958             !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
1959             DO  k = nzb, nzt_diff
1960                !$acc loop vector collapse( 2 ) reduction( +: s1 )
1961                DO  i = nxl, nxr
1962                   DO  j = nys, nyn
1963                      s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )    &
1964                                    * ( ( q(k+1,j,i) - ql(k+1,j,i) ) &
1965                                      - ( q(k,j,i) - ql(k,j,i) ) )   &
1966                                    * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1967                   ENDDO
1968                ENDDO
1969                sums_l(k,51,tn) = s1
1970             ENDDO
1971             !$acc end parallel loop
1972
1973          ENDIF
1974
1975       ENDIF
1976!
1977!--    Passive scalar flux
1978       IF ( passive_scalar )  THEN
1979
1980          !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
1981          DO  k = nzb, nzt_diff
1982             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1983             DO  i = nxl, nxr
1984                DO  j = nys, nyn
1985                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
1986                                 * ( q(k+1,j,i) - q(k,j,i) )     &
1987                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1988                ENDDO
1989             ENDDO
1990             sums_l(k,48,tn) = s1
1991          ENDDO
1992          !$acc end parallel loop
1993
1994       ENDIF
1995
1996       IF ( use_surface_fluxes )  THEN
1997
1998          !$OMP DO
1999          !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
2000          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2001          DO  i = nxl, nxr
2002             DO  j =  nys, nyn
2003!
2004!--             Subgridscale fluxes in the Prandtl layer
2005                s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
2006                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
2007                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
2008                s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
2009                s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
2010             ENDDO
2011          ENDDO
2012          sums_l(nzb,12,tn) = s1
2013          sums_l(nzb,14,tn) = s2
2014          sums_l(nzb,16,tn) = s3
2015          sums_l(nzb,58,tn) = s4
2016          sums_l(nzb,61,tn) = s5
2017          !$acc end parallel
2018
2019          IF ( ocean )  THEN
2020
2021             !$OMP DO
2022             !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
2023             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2024             DO  i = nxl, nxr
2025                DO  j =  nys, nyn
2026                   s1 = s1 + saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
2027                ENDDO
2028             ENDDO
2029             sums_l(nzb,65,tn) = s1
2030             !$acc end parallel
2031
2032          ENDIF
2033
2034          IF ( humidity )  THEN
2035
2036             !$OMP DO
2037             !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
2038             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2039             DO  i = nxl, nxr
2040                DO  j =  nys, nyn
2041                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2042                   s2 = s2 + ( ( 1.0 + 0.61 * q(nzb,j,i) ) * shf(j,i) &
2043                               + 0.61 * pt(nzb,j,i) * qsws(j,i) )
2044                ENDDO
2045             ENDDO
2046             sums_l(nzb,48,tn) = s1
2047             sums_l(nzb,45,tn) = s2
2048             !$acc end parallel
2049
2050             IF ( cloud_droplets )  THEN
2051
2052                !$OMP DO
2053                !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
2054                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2055                DO  i = nxl, nxr
2056                   DO  j =  nys, nyn
2057                      s1 = s1 + ( ( 1.0 + 0.61 * q(nzb,j,i) - ql(nzb,j,i) ) * &
2058                                  shf(j,i) + 0.61 * pt(nzb,j,i) * qsws(j,i) )
2059                   ENDDO
2060                ENDDO
2061                sums_l(nzb,45,tn) = s1
2062                !$acc end parallel
2063
2064             ENDIF
2065
2066             IF ( cloud_physics )  THEN
2067
2068                !$OMP DO
2069                !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2070                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2071                DO  i = nxl, nxr
2072                   DO  j =  nys, nyn
2073!
2074!--                   Formula does not work if ql(nzb) /= 0.0
2075                      s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
2076                   ENDDO
2077                ENDDO
2078                sums_l(nzb,51,tn) = s1
2079                !$acc end parallel
2080
2081             ENDIF
2082
2083          ENDIF
2084
2085          IF ( passive_scalar )  THEN
2086
2087             !$OMP DO
2088             !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2089             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2090             DO  i = nxl, nxr
2091                DO  j =  nys, nyn
2092                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2093                ENDDO
2094             ENDDO
2095             sums_l(nzb,48,tn) = s1
2096             !$acc end parallel
2097
2098          ENDIF
2099
2100       ENDIF
2101
2102!
2103!--    Subgridscale fluxes at the top surface
2104       IF ( use_top_fluxes )  THEN
2105
2106          !$OMP DO
2107          !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
2108          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2109          DO  i = nxl, nxr
2110             DO  j =  nys, nyn
2111                s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
2112                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
2113                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
2114                s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
2115                s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
2116             ENDDO
2117          ENDDO
2118          sums_l(nzt:nzt+1,12,tn) = s1
2119          sums_l(nzt:nzt+1,14,tn) = s2
2120          sums_l(nzt:nzt+1,16,tn) = s3
2121          sums_l(nzt:nzt+1,58,tn) = s4
2122          sums_l(nzt:nzt+1,61,tn) = s5
2123          !$acc end parallel
2124
2125          IF ( ocean )  THEN
2126
2127             !$OMP DO
2128             !$acc parallel present( rmask, saswst, sums_l ) create( s1 )
2129             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2130             DO  i = nxl, nxr
2131                DO  j =  nys, nyn
2132                   s1 = s1 + saswst(j,i) * rmask(j,i,sr)  ! w"sa"
2133                ENDDO
2134             ENDDO
2135             sums_l(nzt,65,tn) = s1
2136             !$acc end parallel
2137
2138          ENDIF
2139
2140          IF ( humidity )  THEN
2141
2142             !$OMP DO
2143             !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
2144             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2145             DO  i = nxl, nxr
2146                DO  j =  nys, nyn
2147                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2148                   s2 = s2 + ( ( 1.0 + 0.61 * q(nzt,j,i) ) * tswst(j,i) + &
2149                               0.61 * pt(nzt,j,i) * qswst(j,i) )
2150                ENDDO
2151             ENDDO
2152             sums_l(nzt,48,tn) = s1
2153             sums_l(nzt,45,tn) = s2
2154             !$acc end parallel
2155
2156             IF ( cloud_droplets )  THEN
2157
2158                !$OMP DO
2159                !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
2160                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2161                DO  i = nxl, nxr
2162                   DO  j =  nys, nyn
2163                      s1 = s1 + ( ( 1.0 + 0.61 * q(nzt,j,i) - ql(nzt,j,i) ) * &
2164                                  tswst(j,i) + 0.61 * pt(nzt,j,i) * qswst(j,i) )
2165                   ENDDO
2166                ENDDO
2167                sums_l(nzt,45,tn) = s1
2168                !$acc end parallel
2169
2170             ENDIF
2171
2172             IF ( cloud_physics )  THEN
2173
2174                !$OMP DO
2175                !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2176                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2177                DO  i = nxl, nxr
2178                   DO  j =  nys, nyn
2179!
2180!--                   Formula does not work if ql(nzb) /= 0.0
2181                      s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2182                   ENDDO
2183                ENDDO
2184                sums_l(nzt,51,tn) = s1
2185                !$acc end parallel
2186
2187             ENDIF
2188
2189          ENDIF
2190
2191          IF ( passive_scalar )  THEN
2192
2193             !$OMP DO
2194             !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2195             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2196             DO  i = nxl, nxr
2197                DO  j =  nys, nyn
2198                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2199                ENDDO
2200             ENDDO
2201             sums_l(nzt,48,tn) = s1
2202             !$acc end parallel
2203
2204          ENDIF
2205
2206       ENDIF
2207
2208!
2209!--    Resolved fluxes (can be computed for all horizontal points)
2210!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
2211!--    ----  speaking the following k-loop would have to be split up and
2212!--          rearranged according to the staggered grid.
2213       !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
2214       DO  k = nzb, nzt_diff
2215          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2216          DO  i = nxl, nxr
2217             DO  j = nys, nyn
2218                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
2219                              u(k+1,j,i) - hom(k+1,1,1,sr) )
2220                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
2221                              v(k+1,j,i) - hom(k+1,1,2,sr) )
2222                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
2223                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
2224!
2225!--             Higher moments
2226                s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2227                s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2228!
2229!--             Energy flux w*e* (has to be adjusted?)
2230                s3 = s3 + w(k,j,i) * 0.5 * ( ust**2 + vst**2 + w(k,j,i)**2 ) &
2231                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2232             ENDDO
2233          ENDDO
2234          sums_l(k,35,tn) = s1
2235          sums_l(k,36,tn) = s2
2236          sums_l(k,37,tn) = s3
2237       ENDDO
2238       !$acc end parallel loop
2239
2240!
2241!--    Salinity flux and density (density does not belong to here,
2242!--    but so far there is no other suitable place to calculate)
2243       IF ( ocean )  THEN
2244
2245          IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
2246
2247             !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
2248             DO  k = nzb, nzt_diff
2249                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2250                DO  i = nxl, nxr
2251                   DO  j = nys, nyn
2252                      s1 = s1 + 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) +   &
2253                                        sa(k+1,j,i) - hom(k+1,1,23,sr) ) &
2254                                    * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2255                   ENDDO
2256                ENDDO
2257                sums_l(k,66,tn) = s1
2258             ENDDO
2259             !$acc end parallel loop
2260
2261          ENDIF
2262
2263          !$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
2264          DO  k = nzb, nzt_diff
2265             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2266             DO  i = nxl, nxr
2267                DO  j = nys, nyn
2268                   s1 = s1 + rho(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2269                   s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2270                ENDDO
2271             ENDDO
2272             sums_l(k,64,tn) = s1
2273             sums_l(k,71,tn) = s2
2274          ENDDO
2275          !$acc end parallel loop
2276
2277       ENDIF
2278
2279!
2280!--    Buoyancy flux, water flux, humidity flux, liquid water
2281!--    content, rain drop concentration and rain water content
2282       IF ( humidity )  THEN
2283
2284          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
2285
2286             !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2287             DO  k = nzb, nzt_diff
2288                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2289                DO  i = nxl, nxr
2290                   DO  j = nys, nyn
2291                      s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
2292                                        vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
2293                                      w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2294                   ENDDO
2295                ENDDO
2296                sums_l(k,46,tn) = s1
2297             ENDDO
2298             !$acc end parallel loop
2299
2300             IF ( .NOT. cloud_droplets )  THEN
2301
2302                !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
2303                DO  k = nzb, nzt_diff
2304                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2305                   DO  i = nxl, nxr
2306                      DO  j = nys, nyn
2307                         s1 = s1 + 0.5 * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
2308                                           ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
2309                                       * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2310                      ENDDO
2311                   ENDDO
2312                   sums_l(k,52,tn) = s1
2313                ENDDO
2314                !$acc end parallel loop
2315
2316                IF ( icloud_scheme == 0  )  THEN
2317
2318                   !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
2319                   DO  k = nzb, nzt_diff
2320                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2321                      DO  i = nxl, nxr
2322                         DO  j = nys, nyn
2323                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2324                            s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2325                         ENDDO
2326                      ENDDO
2327                      sums_l(k,54,tn) = s1
2328                      sums_l(k,75,tn) = s2
2329                   ENDDO
2330                   !$acc end parallel loop
2331
2332                   IF ( precipitation )  THEN
2333
2334                      !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2335                      DO  k = nzb, nzt_diff
2336                         !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2337                         DO  i = nxl, nxr
2338                            DO  j = nys, nyn
2339                               s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2340                               s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2341                               s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2342                            ENDDO
2343                         ENDDO
2344                         sums_l(k,73,tn) = s1
2345                         sums_l(k,74,tn) = s2
2346                         sums_l(k,76,tn) = s3
2347                      ENDDO
2348                      !$acc end parallel loop
2349
2350                   ENDIF
2351
2352                ELSE
2353
2354                   !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2355                   DO  k = nzb, nzt_diff
2356                      !$acc loop vector collapse( 2 ) reduction( +: s1 )
2357                      DO  i = nxl, nxr
2358                         DO  j = nys, nyn
2359                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2360                         ENDDO
2361                      ENDDO
2362                      sums_l(k,54,tn) = s1
2363                   ENDDO
2364                   !$acc end parallel loop
2365
2366                ENDIF
2367
2368             ELSE
2369
2370                !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2371                DO  k = nzb, nzt_diff
2372                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2373                   DO  i = nxl, nxr
2374                      DO  j = nys, nyn
2375                         s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2376                      ENDDO
2377                   ENDDO
2378                   sums_l(k,54,tn) = s1
2379                ENDDO
2380                !$acc end parallel loop
2381
2382             ENDIF
2383
2384          ELSE
2385
2386             IF( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2387
2388                !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2389                DO  k = nzb, nzt_diff
2390                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2391                   DO  i = nxl, nxr
2392                      DO  j = nys, nyn
2393                         s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
2394                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
2395                                       * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2396                      ENDDO
2397                   ENDDO
2398                   sums_l(k,46,tn) = s1
2399                ENDDO
2400                !$acc end parallel loop
2401
2402             ELSEIF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
2403
2404                !$acc parallel loop present( hom, sums_l )
2405                DO  k = nzb, nzt_diff
2406                   sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
2407                                             0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
2408                ENDDO
2409                !$acc end parallel loop
2410
2411             ENDIF
2412
2413          ENDIF
2414
2415       ENDIF
2416!
2417!--    Passive scalar flux
2418       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
2419
2420          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2421          DO  k = nzb, nzt_diff
2422             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2423             DO  i = nxl, nxr
2424                DO  j = nys, nyn
2425                   s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
2426                                     q(k+1,j,i) - hom(k+1,1,41,sr) ) &
2427                                 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2428                ENDDO
2429             ENDDO
2430             sums_l(k,49,tn) = s1
2431          ENDDO
2432          !$acc end parallel loop
2433
2434       ENDIF
2435
2436!
2437!--    For speed optimization fluxes which have been computed in part directly
2438!--    inside the WS advection routines are treated seperatly
2439!--    Momentum fluxes first:
2440       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
2441
2442          !$OMP DO
2443          !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2444          DO  k = nzb, nzt_diff
2445             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2446             DO  i = nxl, nxr
2447                DO  j = nys, nyn
2448                   ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
2449                               u(k+1,j,i) - hom(k+1,1,1,sr) )
2450                   vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
2451                               v(k+1,j,i) - hom(k+1,1,2,sr) )
2452!
2453!--                Momentum flux w*u*
2454                   s1 = s1 + 0.5 * ( w(k,j,i-1) + w(k,j,i) ) &
2455                                 * ust * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2456!
2457!--                Momentum flux w*v*
2458                   s2 = s2 + 0.5 * ( w(k,j-1,i) + w(k,j,i) ) &
2459                                 * vst * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2460                ENDDO
2461             ENDDO
2462             sums_l(k,13,tn) = s1
2463             sums_l(k,15,tn) = s1
2464          ENDDO
2465          !$acc end parallel loop
2466
2467       ENDIF
2468
2469       IF ( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2470
2471          !$OMP DO
2472          !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
2473          DO  k = nzb, nzt_diff
2474             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2475             DO  i = nxl, nxr
2476                DO  j = nys, nyn
2477!
2478!--                Vertical heat flux
2479                   s1 = s1 + 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
2480                                     pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
2481                                 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2482                ENDDO
2483             ENDDO
2484             sums_l(k,17,tn) = s1
2485          ENDDO
2486          !$acc end parallel loop
2487
2488          IF ( humidity )  THEN
2489
2490             !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2491             DO  k = nzb, nzt_diff
2492                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2493                DO  i = nxl, nxr
2494                   DO  j = nys, nyn
2495                      s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
2496                                        q(k+1,j,i) - hom(k+1,1,41,sr) ) &
2497                                    * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2498                   ENDDO
2499                ENDDO
2500                sums_l(k,49,tn) = s1
2501             ENDDO
2502             !$acc end parallel loop
2503
2504          ENDIF
2505
2506       ENDIF
2507
2508
2509!
2510!--    Density at top follows Neumann condition
2511       IF ( ocean )  THEN
2512          !$acc parallel present( sums_l )
2513          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
2514          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
2515          !$acc end parallel
2516       ENDIF
2517
2518!
2519!--    Divergence of vertical flux of resolved scale energy and pressure
2520!--    fluctuations as well as flux of pressure fluctuation itself (68).
2521!--    First calculate the products, then the divergence.
2522!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
2523       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
2524
2525          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
2526          sums_ll = 0.0  ! local array
2527
2528          !$OMP DO
2529          DO  i = nxl, nxr
2530             DO  j = nys, nyn
2531                DO  k = nzb_s_inner(j,i)+1, nzt
2532
2533                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
2534                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
2535                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
2536                           ) )**2                                          &
2537                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
2538                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
2539                           ) )**2                                          &
2540                   + w(k,j,i)**2                                  )
2541
2542                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
2543                                               * ( p(k,j,i) + p(k+1,j,i) )
2544
2545                ENDDO
2546             ENDDO
2547          ENDDO
2548          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
2549          sums_ll(nzt+1,1) = 0.0
2550          sums_ll(0,2)     = 0.0
2551          sums_ll(nzt+1,2) = 0.0
2552
2553          DO  k = nzb+1, nzt
2554             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
2555             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
2556             sums_l(k,68,tn) = sums_ll(k,2)
2557          ENDDO
2558          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
2559          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
2560          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
2561
2562       ENDIF
2563
2564!
2565!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
2566       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
2567
2568          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
2569          !$OMP DO
2570          DO  i = nxl, nxr
2571             DO  j = nys, nyn
2572                DO  k = nzb_s_inner(j,i)+1, nzt
2573
2574                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
2575                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2576                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
2577                                                             ) * ddzw(k)
2578
2579                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
2580                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2581                                                              )
2582
2583                ENDDO
2584             ENDDO
2585          ENDDO
2586          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
2587          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
2588
2589       ENDIF
2590
2591!
2592!--    Horizontal heat fluxes (subgrid, resolved, total).
2593!--    Do it only, if profiles shall be plotted.
2594       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
2595
2596          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
2597          !$OMP DO
2598          DO  i = nxl, nxr
2599             DO  j = nys, nyn
2600                DO  k = nzb_s_inner(j,i)+1, nzt
2601!
2602!--                Subgrid horizontal heat fluxes u"pt", v"pt"
2603                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
2604                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
2605                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
2606                                                 * ddx * rmask(j,i,sr)
2607                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
2608                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
2609                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
2610                                                 * ddy * rmask(j,i,sr)
2611!
2612!--                Resolved horizontal heat fluxes u*pt*, v*pt*
2613                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
2614                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
2615                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
2616                                                 pt(k,j,i)   - hom(k,1,4,sr) )
2617                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
2618                                 pt(k,j,i)   - hom(k,1,4,sr) )
2619                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
2620                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
2621                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
2622                                                 pt(k,j,i)   - hom(k,1,4,sr) )
2623                ENDDO
2624             ENDDO
2625          ENDDO
2626!
2627!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
2628          sums_l(nzb,58,tn) = 0.0
2629          sums_l(nzb,59,tn) = 0.0
2630          sums_l(nzb,60,tn) = 0.0
2631          sums_l(nzb,61,tn) = 0.0
2632          sums_l(nzb,62,tn) = 0.0
2633          sums_l(nzb,63,tn) = 0.0
2634
2635       ENDIF
2636
2637!
2638!--    Calculate the user-defined profiles
2639       CALL user_statistics( 'profiles', sr, tn )
2640       !$OMP END PARALLEL
2641
2642!
2643!--    Summation of thread sums
2644       IF ( threads_per_task > 1 )  THEN
2645          STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
2646          DO  i = 1, threads_per_task-1
2647             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
2648             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
2649             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
2650                                      sums_l(:,45:pr_palm,i)
2651             IF ( max_pr_user > 0 )  THEN
2652                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
2653                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
2654                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
2655             ENDIF
2656          ENDDO
2657       ENDIF
2658
2659       !$acc update host( hom, sums, sums_l )
2660
2661#if defined( __parallel )
2662
2663!
2664!--    Compute total sum from local sums
2665       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2666       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
2667                           MPI_SUM, comm2d, ierr )
2668#else
2669       sums = sums_l(:,:,0)
2670#endif
2671
2672!
2673!--    Final values are obtained by division by the total number of grid points
2674!--    used for summation. After that store profiles.
2675!--    Profiles:
2676       DO  k = nzb, nzt+1
2677          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
2678          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
2679          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
2680          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
2681          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
2682          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
2683          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
2684          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
2685          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
2686          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
2687          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
2688          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
2689          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
2690          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
2691       ENDDO
2692
2693!--    Upstream-parts
2694       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
2695!--    u* and so on
2696!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
2697!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
2698!--    above the topography, they are being divided by ngp_2dh(sr)
2699       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
2700                                    ngp_2dh(sr)
2701       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
2702                                    ngp_2dh(sr)
2703!--    eges, e*
2704       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
2705                                    ngp_3d(sr)
2706!--    Old and new divergence
2707       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
2708                                    ngp_3d_inner(sr)
2709
2710!--    User-defined profiles
2711       IF ( max_pr_user > 0 )  THEN
2712          DO  k = nzb, nzt+1
2713             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
2714                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
2715                                    ngp_2dh_s_inner(k,sr)
2716          ENDDO
2717       ENDIF
2718
2719!
2720!--    Collect horizontal average in hom.
2721!--    Compute deduced averages (e.g. total heat flux)
2722       hom(:,1,3,sr)  = sums(:,3)      ! w
2723       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
2724       hom(:,1,9,sr)  = sums(:,9)      ! km
2725       hom(:,1,10,sr) = sums(:,10)     ! kh
2726       hom(:,1,11,sr) = sums(:,11)     ! l
2727       hom(:,1,12,sr) = sums(:,12)     ! w"u"
2728       hom(:,1,13,sr) = sums(:,13)     ! w*u*
2729       hom(:,1,14,sr) = sums(:,14)     ! w"v"
2730       hom(:,1,15,sr) = sums(:,15)     ! w*v*
2731       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
2732       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
2733       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
2734       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
2735       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
2736       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
2737       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
2738                                       ! profile 24 is initial profile (sa)
2739                                       ! profiles 25-29 left empty for initial
2740                                       ! profiles
2741       hom(:,1,30,sr) = sums(:,30)     ! u*2
2742       hom(:,1,31,sr) = sums(:,31)     ! v*2
2743       hom(:,1,32,sr) = sums(:,32)     ! w*2
2744       hom(:,1,33,sr) = sums(:,33)     ! pt*2
2745       hom(:,1,34,sr) = sums(:,34)     ! e*
2746       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
2747       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
2748       hom(:,1,37,sr) = sums(:,37)     ! w*e*
2749       hom(:,1,38,sr) = sums(:,38)     ! w*3
2750       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
2751       hom(:,1,40,sr) = sums(:,40)     ! p
2752       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
2753       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
2754       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
2755       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
2756       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
2757       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
2758       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
2759       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
2760       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
2761       hom(:,1,54,sr) = sums(:,54)     ! ql
2762       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
2763       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
2764       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
2765       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
2766       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
2767       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
2768       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
2769       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
2770       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
2771       hom(:,1,64,sr) = sums(:,64)     ! rho
2772       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
2773       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
2774       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
2775       hom(:,1,68,sr) = sums(:,68)     ! w*p*
2776       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
2777       hom(:,1,70,sr) = sums(:,70)     ! q*2
2778       hom(:,1,71,sr) = sums(:,71)     ! prho
2779       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
2780       hom(:,1,73,sr) = sums(:,73)     ! nr
2781       hom(:,1,74,sr) = sums(:,74)     ! qr
2782       hom(:,1,75,sr) = sums(:,75)     ! qc
2783       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
2784                                       ! 77 is initial density profile
2785       hom(:,1,78,sr) = ug             ! ug
2786       hom(:,1,79,sr) = vg             ! vg
2787       hom(:,1,80,sr) = w_subs         ! w_subs
2788
2789       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
2790                                       ! upstream-parts u_x, u_y, u_z, v_x,
2791                                       ! v_y, usw. (in last but one profile)
2792       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
2793                                       ! u*, w'u', w'v', t* (in last profile)
2794
2795       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
2796          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
2797                               sums(:,pr_palm+1:pr_palm+max_pr_user)
2798       ENDIF
2799
2800!
2801!--    Determine the boundary layer height using two different schemes.
2802!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
2803!--    first relative minimum (maximum) of the total heat flux.
2804!--    The corresponding height is assumed as the boundary layer height, if it
2805!--    is less than 1.5 times the height where the heat flux becomes negative
2806!--    (positive) for the first time.
2807       z_i(1) = 0.0
2808       first = .TRUE.
2809
2810       IF ( ocean )  THEN
2811          DO  k = nzt, nzb+1, -1
2812             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
2813                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
2814                first = .FALSE.
2815                height = zw(k)
2816             ENDIF
2817             IF ( hom(k,1,18,sr) < 0.0  .AND. &
2818                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
2819                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
2820                IF ( zw(k) < 1.5 * height )  THEN
2821                   z_i(1) = zw(k)
2822                ELSE
2823                   z_i(1) = height
2824                ENDIF
2825                EXIT
2826             ENDIF
2827          ENDDO
2828       ELSE
2829          DO  k = nzb, nzt-1
2830             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
2831                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
2832                first = .FALSE.
2833                height = zw(k)
2834             ENDIF
2835             IF ( hom(k,1,18,sr) < 0.0  .AND. &
2836                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
2837                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
2838                IF ( zw(k) < 1.5 * height )  THEN
2839                   z_i(1) = zw(k)
2840                ELSE
2841                   z_i(1) = height
2842                ENDIF
2843                EXIT
2844             ENDIF
2845          ENDDO
2846       ENDIF
2847
2848!
2849!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
2850!--    by Uhlenbrock(2006). The boundary layer height is the height with the
2851!--    maximal local temperature gradient: starting from the second (the last
2852!--    but one) vertical gridpoint, the local gradient must be at least
2853!--    0.2K/100m and greater than the next four gradients.
2854!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
2855!--             ocean case!
2856       z_i(2) = 0.0
2857       DO  k = nzb+1, nzt+1
2858          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
2859       ENDDO
2860       dptdz_threshold = 0.2 / 100.0
2861
2862       IF ( ocean )  THEN
2863          DO  k = nzt+1, nzb+5, -1
2864             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
2865                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
2866                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
2867                z_i(2) = zw(k-1)
2868                EXIT
2869             ENDIF
2870          ENDDO
2871       ELSE
2872          DO  k = nzb+1, nzt-3
2873             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
2874                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
2875                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
2876                z_i(2) = zw(k-1)
2877                EXIT
2878             ENDIF
2879          ENDDO
2880       ENDIF
2881
2882       hom(nzb+6,1,pr_palm,sr) = z_i(1)
2883       hom(nzb+7,1,pr_palm,sr) = z_i(2)
2884
2885!
2886!--    Computation of both the characteristic vertical velocity and
2887!--    the characteristic convective boundary layer temperature.
2888!--    The horizontal average at nzb+1 is input for the average temperature.
2889       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
2890           .AND.  z_i(1) /= 0.0 )  THEN
2891          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
2892                                       hom(nzb,1,18,sr) *      &
2893                                       ABS( z_i(1) ) )**0.333333333
2894!--       so far this only works if Prandtl layer is used
2895          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
2896       ELSE
2897          hom(nzb+8,1,pr_palm,sr)  = 0.0
2898          hom(nzb+11,1,pr_palm,sr) = 0.0
2899       ENDIF
2900
2901!
2902!--    Collect the time series quantities
2903       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
2904       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
2905       ts_value(3,sr) = dt_3d
2906       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
2907       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
2908       ts_value(6,sr) = u_max
2909       ts_value(7,sr) = v_max
2910       ts_value(8,sr) = w_max
2911       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
2912       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
2913       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
2914       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
2915       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
2916       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
2917       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
2918       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
2919       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
2920       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
2921       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
2922       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
2923       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
2924
2925       IF ( ts_value(5,sr) /= 0.0 )  THEN
2926          ts_value(22,sr) = ts_value(4,sr)**2 / &
2927                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
2928       ELSE
2929          ts_value(22,sr) = 10000.0
2930       ENDIF
2931
2932       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
2933!
2934!--    Calculate additional statistics provided by the user interface
2935       CALL user_statistics( 'time_series', sr, 0 )
2936
2937    ENDDO    ! loop of the subregions
2938
2939    !$acc end data
2940
2941!
2942!-- If required, sum up horizontal averages for subsequent time averaging
2943    IF ( do_sum )  THEN
2944       IF ( average_count_pr == 0 )  hom_sum = 0.0
2945       hom_sum = hom_sum + hom(:,1,:,:)
2946       average_count_pr = average_count_pr + 1
2947       do_sum = .FALSE.
2948    ENDIF
2949
2950!
2951!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
2952!-- This flag is reset after each time step in time_integration.
2953    flow_statistics_called = .TRUE.
2954
2955    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
2956
2957
2958 END SUBROUTINE flow_statistics
2959#endif
Note: See TracBrowser for help on using the repository browser.