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

Last change on this file since 1555 was 1555, checked in by maronga, 6 years ago

LSM output of r_a and r_s added

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