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

Last change on this file since 1499 was 1499, checked in by suehring, 7 years ago

last commit documented

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