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

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

last commit documented

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