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

Last change on this file since 1255 was 1242, checked in by heinze, 11 years ago

Last commmit documented

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