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

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

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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