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

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

New:


openACC porting of reduction operations
additional 3D-flag arrays for replacing the 2D-index arrays nzb_s_inner and nzb_diff_s_inner
(flow_statistics, init_grid, init_3d_model, modules, palm, pres, time_integration)

Changed:


for PGI/openACC performance reasons set default compile options for openACC to "-ta=nocache",
and set environment variable PGI_ACC_SYNCHRONOUS=1
(MAKE.inc.pgi.openacc, palm_simple_run)

wall_flags_0 changed to 32bit INTEGER, additional array wall_flags_00 introduced to hold
bits 32-63
(advec_ws, init_grid, modules, palm)

Errors:


dummy argument tri in 1d-routines replaced by tri_for_1d because of name
conflict with arry tri in module arrays_3d
(tridia_solver)

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