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

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

land surface model released

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