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

Last change on this file since 1222 was 1222, checked in by raasch, 8 years ago

last commit documented

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