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

Last change on this file since 1374 was 1374, checked in by raasch, 7 years ago

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

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