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

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

minor bugfixes in lasm/radiation

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