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

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

last commit documented

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