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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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