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

Last change on this file since 1239 was 1239, checked in by heinze, 8 years ago

routines for nudging and large scale forcing from external file added

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