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

Last change on this file since 1319 was 1319, checked in by raasch, 7 years ago

last commit documented

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