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

Last change on this file since 1567 was 1567, checked in by suehring, 6 years ago

Bugfixes in monotonic limter.

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