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

Last change on this file since 1383 was 1383, checked in by boeske, 7 years ago

last commit documented

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